Haploid assembly method based on deep learning

Abstract

The invention discloses a haplotype assembly method based on deep learning. According to the method, the readings are clustered by learning the correlation among the sequencing readings, so that haplotype assembly is realized. Feature representation of learning readings is based on a new large language model autoregressive infrastructure RetNet, and dependence between readings which are far away from a reference genome comparison position can be learned through a multi-scale retention mechanism of the new large language model autoregressive infrastructure RetNet. Based on reading feature representation learned by a RetNet model, a SpectraNet model is used to realize a clustering process of sequencing reading, and finally a haplotype is constructed according to a reading cluster. The method is simple and easy to use, in experiments of simulation data and real data, the method is obviously good in performance in diploid and polyploid haplotype assembly tasks based on long readings or short readings, and compared with other haplotype assembly tools, the method has higher accuracy and continuity.

Description

Technical Field

[0001] The present invention relates to the field of haplotype assembly in bioinformatics, and in particular to a haplotype assembly method based on deep learning. Background Art

[0002] The human genome consists of two sets of homologous chromosomes, which come from the father and the mother respectively. There are a certain number of sequence variations between the two. The combination of different genetic loci has an important impact on the biological phenotype. Haplotype refers to a combination of a group of related single nucleotide polymorphism (SNP) alleles located on a chromosome or in a certain region. Diploid individuals have two haplotypes, and two sets of sequence sets can be obtained through haplotype assembly. Genome assembly refers to the process of using DNA fragments (i.e. reads) obtained by sequencing technology to splice these fragments into longer and continuous sequences (such as contigs and scaffolds) through computational methods and bioinformatics tools, and finally construct a nearly complete or complete genome sequence. Haplotype assembly focuses more on further analyzing the genetic variation and gene combination on a single chromatid (i.e. haploid) in an organism based on genome assembly. Since organisms are usually genetically diploid or polyploid, that is, there are two or more alleles at each gene locus, the goal of haplotype assembly is to reconstruct the complete or nearly complete haplotype sequence on each chromatid of the organism based on the reads obtained by sequencing and known genetic variation information through computational methods and bioinformatics algorithms.

[0003] The application of haplotype assembly is very useful in many research fields. First, haplotype assembly helps to analyze differences between alleles, track individual kinship, and understand biological migration patterns and evolutionary history in population genetics research. Second, haplotype assembly helps to discover excellent allele variations and explore the theory of hybrid vigor in the agricultural field. It plays a very important role in genetic breeding of major crops such as rice, corn, potatoes, and wheat, as well as other plants such as strawberries and litchi. Third, haplotype assembly helps to explore pathogenic mechanisms, dig out pathogenic genes, and find new methods for treating diseases in the research of biological theory and medicine. Researchers used haplotype technology to analyze the causes of skin acne, cerebral palsy, and deafness in patients, and finally found that they were all autosomal recessive genetic diseases caused by heterozygous variations on chromosomes. In addition, haplotype genome sequencing of the fetus can be used to detect whether it has potential genetic diseases. In view of the above three problems that haplotypes can solve, it can be said that obtaining a high-quality haplotype genome of a species will promote the more profound development of industries related to the species, and will have epoch-making significance for promoting plant molecular breeding, accelerating the selection of new varieties, and preventing, diagnosing and researching human diseases.

[0004] DNA sequencing technology is widely used, such as genomics, transcriptomics, epigenomics, etc. The progress of sequencing technology has also promoted the development of haplotype assembly problems, but the limitations of sequencing platforms are still challenging for haplotype assembly problems. At present, sequencing technology can be divided into second-generation sequencing technology (also known as high-throughput sequencing technology) and third-generation sequencing technology according to its sequencing read length. At present, the mainstream second-generation sequencing platform mainly comes from companies such as Illumina and BGI Genomics. It has lower cost, faster sequencing speed and higher accuracy. Due to the short reading length, the number of covered variant sites is small, which is not enough to provide sufficient evidence to completely phase haplotypes, resulting in low continuity of haplotype assembly. The third-generation sequencing technology is represented by single molecule real-time sequencing technology (SMRT, Single Molecule Real-Time) and nanopore sequencing technology developed by PacBio and Oxford Nanopore respectively. The third-generation technology has higher sequencing speed and accuracy. SMRT sequencing can produce long read sequences with an average length of 10-15kb. Although the error rate is relatively high, about 15%, it can be effectively corrected through multiple sequencing. The average read length of Nanopore sequencing technology ranges from 20-50kb, and the longest can reach the Mb level. The sequencing accuracy is similar to that of PacBio, which is about 86%, and the errors are mainly random sequencing errors. If both the positive and negative strands of DNA are sequenced, the accuracy can reach about 96%. PacBio also launched HiFi sequencing in 2019, which can provide base-level resolution and 99.9% single-molecule reading accuracy, which is comparable to short-read sequencing and Sanger sequencing. The third-generation sequencing technology has longer read lengths and higher accuracy, so it can provide better genome coverage and lower error rates, which can greatly improve the quality of genome assembly.

[0005] The input data for calculating the haplotype assembly problem is a collection of a large number of reads, each of which comes from an unknown haplotype. Haplotype assembly strategies using sequencing reads as input are mainly divided into two categories: haplotype assembly based on a reference genome and de novo haplotype assembly.

[0006] Reference-based haplotype assembly methods assume that reference sequences, aligned reads, and called variants can be used as inputs, and all reads are divided into k subsets (k is ploidy), and finally the corresponding k different haplotypes are assembled. This type of method mainly assembles based on grouped reads, but factors such as sequencing errors, read length, and coverage of reads make haplotype assembly challenging, and with the increase of biological ploidy, the difficulty of accurately assembling haplotypes increases. As mentioned above, in order to achieve k-division of reads and assemble accurate haplotypes, the noisy parts must be corrected. To this end, researchers have proposed several multiple combination optimization models. These include minimum fragment removal (MFR), minimum SNP removal (MSR), minimum fragment cutting (MFC), and minimum error correction (MEC). Since most of these models have been proven to be NP-hard, researchers have proposed applying heuristic strategies for haplotype assembly. Representative methods include HapCUT, HapCUT2, and Whatshap, etc. However, these methods can only solve the problem of diploid haplotype assembly. Methods that can handle haplotype assembly of diploid and polyploid genomes include: HapTree, Ranbow, and H-PoP. In recent years, there are several methods for haplotype assembly of polyploid species or viral quasispecies using deep learning models. GAEseq and CAECseq are two haplotype assembly models based on graph autoencoders and convolutional autoencoders, respectively. XHap uses the attention mechanism of Transformer to learn the correlation between reads. However, the Kernel k-means clustering algorithm is used to simply cluster the reads based on the correlation between the learned features.

[0007] Many tools have been developed for de novo haplotype assembly methods. For example, Dipasm uses only two types of input data (HiFi and Hi-C) and does not rely on parental sequences. It can generate haploid assembled genomes without using reference sequences, but it is prone to incorrect division of highly heterozygous regions. Hifiasm is an assembly tool based on the OLC algorithm that supports assembly using single or multiple long-read data sets. Phasebook generates high-quality haplotype-aware genome assemblies based on long-read data. WHdenovo is a de novo assembly method that uses pedigree sequence graphs to perform haplotype-aware de novo assembly of related diploid individuals. ALLHiC is a new auxiliary assembly software developed for polyploid species, and it uses it to complete the haploid genome assembly of tetraploid sugarcane based on heterozygous site information. FALCON-Unzip is the first method to use long-read sequencing technology (PacBio) for de novo assembly of haploid sequences. FALCON-Unzip is based on the initial assembly of FALCON, by analyzing the haplotype information on PacBio reads, looking for heterozygosity differences, such as SNPs in Contigs, and then typing the heterozygosity difference sequences on the Primary Contig, reassembling them into haplotype sequence (haplotig) blocks, and finally obtaining a diploid genome assembly. HiCanu is an improved version based on Canu, mainly optimized for PacBioHiFi data.

[0008] At present, although haplotype assembly based on reference genome and de novo haplotype assembly methods have achieved good results, there are still several unresolved issues that need further study:

[0009] (1) Traditional haplotype assembly methods based on reference genomes, such as heuristic algorithms, cluster reads based on the overlap between reads aligned to the reference genome and the conflict and consistency of SNP information covered by the reads. However, the clustering results are affected by sequencing technology. Short-read sequencing has higher accuracy, but the read length is shorter, resulting in a limited number of SNPs covered by each read. Long-read sequencing can provide longer reads at a lower cost, but the error rate is higher.

[0010] (2) Factors such as sequencing read coverage and biological ploidy also affect haplotype assembly results. For low-coverage genomic regions, assembly accuracy is reduced. Although some existing haplotype assembly methods have high assembly accuracy, they are only applicable to diploid species and have certain limitations in their ability to generalize species ploidy.

[0011] (3) The de novo haplotype assembly method uses supplementary data to phase haplotypes, so it is generally not affected by species heterozygosity and can handle the haplotype assembly problem of polyploid high-heterozygous species. However, the cost of obtaining supplementary data in real experimental scenarios is relatively high.

[0012] The existence of these problems limits the existing haplotype assembly methods from achieving more satisfactory results. Summary of the invention

[0013] The technical problem to be solved by the present invention is to provide a haplotype assembly method based on deep learning, which can take long reads or short reads as input, in view of the deficiencies of the above-mentioned prior art. The method is simple to use and has high continuity and accuracy.

[0014] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0015] A haplotype assembly method based on deep learning, comprising the following steps:

[0016] S1 data preprocessing;

[0017] S11 uses BWA-MEM to align the reads with the reference genome to obtain a SAM file based on the input read data and the reference genome, extracts SNP positions and genotype information based on the SAM file, and constructs a SNPs matrix;

[0018] S12 Since sequencing reads covering only one SNP site do not provide any information required for haplotype phasing, sequencing reads covering only one SNP site are removed from the constructed SNPs matrix;

[0019] S13 performs One-Hot encoding on the SNPs matrix, converting the one-dimensional sequencing read allele sequence into two-dimensional read data t i ;

[0020] S2 obtains reading embeddings;

[0021] One-Hot encoded sequencing reads are transformed into i Map to the specified lower dimension to obtain the low-dimensional potential representation h i ;

[0022] S3 extracts reading features;

[0023] S31 will learn the potential representation h i Stacked together to form a reading embedding matrix, filling a virtual dimension to form a three-dimensional tensor T 0 ;

[0024] S32 T 0Input to the RetNet model for feature extraction to learn the feature representation of sequencing reads, and design constraints based on the loss function to indirectly optimize the MEC score;

[0025] S4 Read clustering;

[0026] The SpectralNet model is used as the clustering layer to learn the similarity between sequencing read pairs based on the feature representation extracted by the RetNet model to complete the clustering process of sequencing reads;

[0027] S5 cluster optimization and obtain the resulting haplotype;

[0028] S51 aims to optimize the loss function and obtains the reading clustering result by iteratively running the previous steps;

[0029] S52 further refines the clustering results using an optimization algorithm. Finally, a haplotype consensus is obtained by majority voting based on all reads in the same cluster. Each haplotype consensus result refers to a haplotype.

[0030] The step S13 is specifically as follows:

[0031] First, the alleles in the SNPs matrix M are one-hot encoded. ij ∈{A,T,C,G}, so the four possible alleles contained in the reads in the matrix M are mapped one by one to the One-Hot encoding, as follows: A→{1,0,0,0}, C→{0,1,0,0}, T→{0,0,0,1}, T→{0,0,1,0}. If a SNP site is not covered by the sequencing reads, the missing allele in the matrix M is encoded as {0,0,0,0}.

[0032] The step S2 is specifically as follows:

[0033] Using the convolutional encoder from CAECSeq, the one-hot encoded sequencing reads t i Map to the specified lower dimension to obtain the low-dimensional potential representation h i . The convolutional encoder contains three two-dimensional convolutional layers and one fully connected layer. The convolutional layer can capture local features by applying a set of learnable convolution kernels (also called convolution kernels or weights) to the input data. By operating on the local area of ​​the input data, the spatial hierarchical structure of the data is retained, and the output feature map can characterize the characteristics of each part of the input, which is conducive to analyzing the two-dimensional spatial features of the sequencing read data. The stacking of three convolutional layers allows the model to learn feature hierarchies from low-level to high-level. The feature map is converted into a one-dimensional vector through the Flatten operation and input into the fully connected layer to obtain a low-dimensional potential representation h. i .

[0034] The formulation of each layer of the encoder part is as follows:

[0035] T 0 =t i

[0036]

[0037] Where '*' represents the convolution operation, and σ uses PReLU as the activation function. and Represents the weight and bias of the lth convolutional layer, W1 dense and B1 dense represents the weight and bias of the last fully connected layer, h i is the 128-dimensional potential representation obtained. The sizes of the convolution kernels in each layer are (4, 5), (1, 5), (1, 3), the filter sizes are 32, 64, 128, and the step size is 1. The decoder part is formulated as follows:

[0038]

[0039] is the reconstructed data obtained by the decoder. The convolutional layer design of the decoder is symmetrical with that of the encoder.

[0040] The step S32 is specifically as follows:

[0041] The RetNet model is used for feature extraction, which can effectively learn the potential features of the readings. RetNet is stacked by L identical blocks, and each RetNet block contains two modules: a multi-scale retention (MSR) module and a feed-forward network (FFN) module. The core of multi-scale retention is the retention mechanism. Through multi-scale gated Rentention, each attention head can model multi-scale information and is not affected by each other. Compared with the standard self-attention mechanism, the retention mechanism has two major characteristics: 1. The introduction of position-related exponential decay terms to replace softmax simplifies the calculation, while retaining the information of the previous step in a decayed form. 2. The introduction of complex space to express position information replaces absolute or relative position encoding, which is easy to convert into a cyclic form.

[0042] First, the learned potential representation h i Stacked together to form a reading embedding matrix, filling a virtual dimension to form a three-dimensional tensor X 0 Input into RetNet. The L-layer model is constructed by stacking multi-scale retention modules and feedforward network modules. This patent uses X 0 As input, the L-layer RetNet model calculates the output X 1 :

[0043] Y l =MSR(LN(X l ))+X l

[0044] X l+1 =FFN(LN(Y l ))+Y l

[0045] Where LN(·) is LayerNorm, FFN(X)=gelu(XW1)W2, W1 and W2 are parameter matrices.

[0046] In this model, a parallel Retention mechanism is used. Given an input matrix X, the parallel representation of the Retention layer is defined as follows:

[0047]

[0048] Retention(X)=(QK T ⊙D)V

[0049] is the complex conjugate of Θ, combining causal masking and exponential decay along relative distances into a single matrix D. n , K n is the matrix that needs to be learned.

[0050] The Retention submodule of each layer of RetNet is also divided into h heads, each head uses a different W q , W k , W v Parameters, and each head uses a different γ constant. In our model, h = 4. Then for input X, the output of the MSR layer is:

[0051] γ=1-2 -5-arange(0，h)

[0052] head i =Retention(X,γ i )

[0053] Y=GroupNorm h (ConCat(head1, ..., head h ))

[0054] MSR(X)=(swish(XW G )⊙Y)W O

[0055] Among them, swish is the activation function used to generate the gating threshold. Since each head uses a different γ parameter, the output of each head needs to be normalized and then concat.

[0056] The step S4 is specifically as follows:

[0057] Based on the feature representation of sequencing reads obtained by the RetNet model, this model uses SpectralNet as the clustering module to cluster the reads into k clusters, and then constructs haplotypes by calculating the consensus sequence based on majority voting.

[0058] The entire SpectralNet training can be divided into three parts: 1. Taking features as input, a twin network is used to calculate the affinity matrix using a given distance metric; 2. By enforcing orthogonality while optimizing the clustering target F θ , mapping the input data points into the eigenspace of their associated graph Laplacian matrix; 3. In the embedding space, the k-means algorithm is used to learn the cluster assignments.

[0059] X 3 is the feature of the sequencing reads output by the last layer of RetNet model, and the twin network converts X 3 Each data point x i Mapped into a space Typically, the network is trained to minimize contrast loss. After the twin network training is completed, we use it to construct an affinity matrix A as the basis for training SpectralNet. The specific calculation formula is as follows.

[0060]

[0061] Where c is the marginal value (usually set to 1). σ is the standard deviation of the Gaussian distribution.

[0062] SpectralNet uses a general neural network to implement the mapping F., that is, yi = F. (xi), and its last layer has an orthogonality constraint. The model is trained by minimizing the loss LSpectralNet, which is defined as follows.

[0063]

[0064] Finally, for y1, ..., y m Use k-means to perform clustering and obtain clustering results.

[0065] The step S51 is specifically as follows:

[0066] This patent uses an unsupervised deep learning to extract similarities between reads and then cluster them for haplotype assembly. This is achieved through an iterative training process of the neural network. Specifically, each training epoch includes the following two steps:

[0067] (1) In the previous step, the feature representation and cluster label of each reading can be obtained. Then, the loss function is used to train the convolutional encoder model and RetNet;

[0068] (2) After obtaining the optimized features, SpectralNet is used again to obtain the clustering labels of the readings.

[0069] These two steps were repeated 2000 times to obtain the clustering results.

[0070] The cluster assignments and read clusters obtained from model training are used as input to the optimization algorithm, which outputs the final optimized cluster reads. Finally, the consensus sequence is extracted as the result haplotype by majority voting.

[0071] The overall loss function is L = L con +β2L reg +λL AE , β2 and λ are weight parameters. The main purpose is to encourage the model to learn the characteristics of the reads more effectively so as to accurately cluster the reads from the same haplotype.

[0072] The contrastive loss is used to effectively handle the relationship between readings. The specific form is as follows:

[0073]

[0074] RetNet model output feature X 3 , the inner product matrix is ​​calculated as the similarity matrix S, where S = X 3 ·(X 3 ) T , S ij represents the similarity between the i-th read and the j-th read. The expected similarity matrix E is constructed based on the pseudo-labels obtained by cluster assignment. If reads i and j belong to the same haplotype, then E ij = 1 is considered as a positive sample pair, otherwise E ij = 0 is considered as a negative sample pair. β1 is a preset boundary value hyperparameter used to distinguish the similarity boundary between positive and negative pairs.

[0075] The loss function is constructed to encourage consistency between the similarity matrix output by the model and the expected similarity matrix, thereby indirectly minimizing the MEC score. When read i and read j are from the same haplotype, E ij =1, L con=(1-S ij ) 2 , the goal is to maximize S ij , when S ij When it is close to 1, 1-S ij The expected similarity matrix E and the correlation matrix S calculated by the model are constrained to be the same for two reads from the same haplotype. When read i and read j come from different haplotypes, if the expected similarity E ij =0, L con =max(0,s i -β1) 2 When S ij When it is less than β1, the loss is 0. ij When it exceeds β1, the loss is (S ij -β1) 2 , encourage S ij Keep it below β1. The synchronization constraint expects the similarity matrix E and the similarity matrix S obtained by the model to behave the same for two reads from different haplotypes. The above two constraints can encourage the model to effectively learn the correlation between read pairs.

[0076] The design of the regularized loss function is based on the improvement of XHap's loss function. Since read pairs with overlapping sites can provide more predictions about the correlation of the read pairs, strong constraints should be imposed on read pairs with overlapping sites. First, a weight based on the total number of overlapping sites is assigned to each pair of reads. When there is no overlap between two read pairs, the weight is 0, avoiding the assignment of unreasonable weights to non-overlapping read pairs. Construct an m*m dynamic weight matrix W dy , specifically defined in the following form:

[0077] W dv (i, j) = overlap ij max(len i ,len j ) where overlap i，i Represents the total number of overlapping sites between read i and read j. i Represents the total number of effective sites of read i. This weight matrix is ​​used to emphasize the correlation between two read pairs, because read pairs with overlap can provide stronger evidence support for haplotype phasing.

[0078] Whether a specific sample pair is strongly positively correlated or strongly negatively correlated is reflected by calculating the correlation value of each pair of overlapping readings i and j: sim is the total number of overlapping sites where the two reads have consistent genotypes, k dissim is the total number of overlapping sites with conflicting genotypes, then the correlation between reads i and j is calculated as:

[0079]

[0080] For overlapping read pairs, if read i and read j are from the same haplotype, then C ij A positive value indicates a positive sample pair. If read i and read j come from different haplotypes, then C ij A negative value indicates a negative sample pair. If the two reads i and j have no overlapping sites, then C ij The value of C is 0. ij The calculation of the reading pair source is consistent with the similarity matrix S predicted by the model. Therefore, the regularization loss is defined as follows:

[0081]

[0082] where F ij is the value of the mask matrix at position (i, j).

[0083] The loss function of the convolutional encoder training process is described as follows:

[0084]

[0085] This loss function represents the reconstruction of the data and the mean square error between the original data ti, where m is the number of samples in the batch.

[0086] The step S52 is specifically as follows:

[0087] Based on the obtained reading clusters and the reading labels obtained by clustering, a heuristic algorithm is used to make possible adjustments to the labels of each reading. Specifically, if a reading is assigned a different label than its current one, which can obtain a better MEC score, the label will be changed. The heuristic algorithm repeatedly tries possible label adjustments for each reading. Heuristic algorithms are often prone to falling into local optimal solutions, especially when the solution space is complex or the quality of solutions varies greatly. Through its "annealing" process, simulated annealing has a higher probability of accepting inferior solutions in the early stages of the algorithm. As the "temperature" gradually decreases, the possibility of accepting inferior solutions is gradually reduced, which helps the algorithm jump out of the local optimal solution and explore a wider solution space. Through the global search capability of simulated annealing, the robustness of the heuristic algorithm to initial values ​​or random fluctuations can be improved, so that the algorithm can maintain relatively stable performance on different runs and different input data to find the best MEC score.

[0088] Compared with the prior art, the present invention has the following beneficial effects:

[0089] The present invention discloses a deep learning model based on unsupervised learning to solve the haplotype assembly problem. The large language model autoregressive infrastructure RetNet is used for sequence modeling, and the feature representation of the readings and the correlation of the global readings are learned through Multi-Scale Retention. With the reading features as input, a twin network is used to calculate the reading affinity matrix using a given distance metric, and the input reading data points are mapped to the eigenspace of the Laplacian matrix of its related graph by enforcing orthogonality while optimizing the clustering target. In the embedded space, the k-means algorithm is used to learn the clustering assignment. In short, the deep spectral clustering method is used to directly perform clustering according to the feature representation of the readings, and after further optimizing the clustering, the consistent sequence is extracted as the haplotype according to the readings in each cluster.

[0090] It performs well in experiments with simulated data and real data. In terms of sequencing coverage, it shows strong robustness, whether it is the simulated short read 10×~30× coverage or the long read 80×~100× coverage. In addition, the model can accept short read data, long read data, or even mixed long and short data as input to the model. It can not only solve the haplotype assembly problem of diploids, but also the haplotype assembly problem of polyploids. In experiments with polyploid data, the experimental results are outstanding compared with other tools.

[0091] The present invention is simple and easy to use, and exhibits strong robustness in terms of model performance for different influencing factors, such as sequencing read length, coverage, species ploidy, etc., and has higher accuracy and continuity than other haplotype assembly methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0092] Figure 1 is a flow chart of the method of the present invention; DETAILED DESCRIPTION

[0093] like Figure 1 As shown, the specific implementation process of the present invention is as follows:

[0094] S1 Data Preprocessing

[0095] S11 Data Processing

[0096] This patent uses read data and reference genome as input, uses BWA-MEM to align the reads with the reference genome to obtain a SAM file, extracts SNP positions and genotype information based on the SAM file, and constructs a SNPs matrix.

[0097] Let M be an m*n SNPs matrix, m is the number of sequencing reads, n is the number of SNP sites covered by all sequencing reads, M ij∈{A,T,C,G}. Since sequencing reads covering only one SNP site do not provide any information required for haplotype phasing, we remove sequencing reads covering only one SNP site in the constructed SNPs matrix M.

[0098] S12 Coding

[0099] One-hot encodes the alleles in the matrix M. ij ∈{A,T,C,G}, so the four possible alleles contained in the reads in the matrix M are mapped one by one to the One-Hot encoding, as follows: A→{1,0,0,0}, C→{0,l,0,0}, T→{0,0,0,1}, T→{0,0,l,0}. If a SNP site is not covered by the sequencing reads, the missing allele in the matrix M is encoded as {0,0,0,0}.

[0100] S2 embeds the reads

[0101] Using a convolutional encoder, the one-hot encoded sequencing reads x i Map to the specified lower dimension to obtain the low-dimensional potential representation h i . The convolutional encoder contains three two-dimensional convolutional layers and one fully connected layer. The convolutional layer can capture local features by applying a set of learnable convolution kernels (also called convolution kernels or weights) to the input data. By operating on the local area of ​​the input data, the spatial hierarchical structure of the data is retained, and the output feature map can characterize the characteristics of each part of the input, which is conducive to analyzing the two-dimensional spatial features of the sequencing read data. The stacking of three convolutional layers allows the model to learn feature hierarchies from low-level to high-level. The feature map is converted into a one-dimensional vector through the Flatten operation and input into the fully connected layer to obtain a low-dimensional potential representation h. i .

[0102] The formulation of each layer of the encoder part is as follows:

[0103] T 0 =t i

[0104]

[0105]

[0106] Where '*' represents the convolution operation, and σ uses PReLU as the activation function. and Represents the weight and bias of the lth convolutional layer, W1 dense and B1 dense represents the weight and bias of the last fully connected layer, h iis the 128-dimensional potential representation obtained. The sizes of the convolution kernels in each layer are (4, 5), (1, 5), (1, 3), the filter sizes are 32, 64, 128, and the step size is 1. The decoder part is formulated as follows:

[0107]

[0108] is the reconstructed data obtained by the decoder. The convolutional layer design of the decoder is symmetrical with that of the encoder.

[0109] S3 extracts read features

[0110] The RetNet model is used for feature extraction, which can effectively learn the potential features of the readings. RetNet is stacked by L identical blocks, and each RetNet block contains two modules: a multi-scale retention (MSR) module and a feed-forward network (FFN) module. The core of multi-scale retention is the retention mechanism. Through multi-scale gated Rentention, each attention head can model multi-scale information and is not affected by each other. Compared with the standard self-attention mechanism, the retention mechanism has two major characteristics: 1. The introduction of position-related exponential decay terms to replace softmax simplifies the calculation, while retaining the information of the previous step in a decayed form. 2. The introduction of complex space to express position information replaces absolute or relative position encoding, which is easy to convert into a cyclic form.

[0111] S31 first transforms the learned potential representation h i Stacked together to form a reading embedding matrix, filling a virtual dimension to form a three-dimensional tensor X 0 Input into RetNet.

[0112] S32 This patent uses X 0 As input, the L-layer RetNet model calculates the output X 1 :

[0113] Y l =MSR(LN(X l ))+X l

[0114] X l+1 =FFN(LN(Y l ))+Y l

[0115] Where LN(·) is LayerNorm, FFN(X)=gelu(XW1)W2, W1 and W2 are parameter matrices.

[0116] In this model, a parallel Retention mechanism is used. Given an input matrix X, the parallel representation of the Retention layer is defined as follows:

[0117]

[0118] Retention(X)=(QK T ⊙D)V

[0119] is the complex conjugate of Θ, combining causal masking and exponential decay along relative distances into a single matrix D. n , K n is the matrix that needs to be learned.

[0120] The Retention submodule of each layer of RetNet is also divided into h heads, each head uses a different W q , W k , W v Parameters, and each head uses a different γ constant. In our model, h = 4. Then for input X, the output of the MSR layer is:

[0121] γ=1-2 -5-arange(0，h)

[0122] head i =Retention(X,γ i )

[0123] Y=GroupNorm h (Concat(head1, ..., head h ))

[0124] MSR(X)=(swish(XW G )⊙Y)W O

[0125] Among them, swish is the activation function used to generate the gating threshold. Since each head uses a different γ parameter, the output of each head needs to be normalized and then concat.

[0126] S4 Clustering of reads

[0127] Based on the feature representation of sequencing reads obtained by the RetNet model, this model uses SpectralNet as the clustering module to cluster the reads into k clusters, and then constructs haplotypes by calculating the consensus sequence based on majority voting.

[0128] The entire SpectralNet training can be divided into three parts: 1. Taking features as input, a twin network is used to calculate the affinity matrix using a given distance metric; 2. By enforcing orthogonality while optimizing the clustering target F θ , mapping the input data points into the eigenspace of their associated graph Laplacian matrix; 3. In the embedding space, the k-means algorithm is used to learn the cluster assignments.

[0129] X 3 is the feature of the sequencing reads output by the last layer of RetNet model, and the twin network converts X 3 Each data point x i Mapped into a space Usually, the network is trained to minimize contrast loss. After the twin network training is completed, we use it to construct an affinity matrix A as the basis for training SpectralNet. The specific calculation formula is as follows.

[0130]

[0131] Where c is the marginal value (usually set to 1). σ is the standard deviation of the Gaussian distribution.

[0132] SpectralNet uses a general neural network to achieve the mapping F θ , that is, y i =F θ (x i ), whose last layer has an orthogonality constraint. The model is trained by minimizing the loss LSpectralNet, which is defined as follows.

[0133]

[0134] Finally, for y1, ..., y m Use k-means to perform clustering and obtain clustering results.

[0135] S5 obtained haplotype

[0136] S51 iterates to obtain clustering results and optimizes

[0137] This patent uses an unsupervised deep learning to extract similarities between reads and then cluster them for haplotype assembly. This is achieved through an iterative training process of the neural network. Specifically, each training epoch includes the following two steps:

[0138] (1) In the previous step, the feature representation and cluster label of each reading can be obtained. Then, the loss function is used to train the convolutional encoder model and RetNet;

[0139] (2) After obtaining the optimized features, SpectralNet is used again to obtain the clustering labels of the readings.

[0140] These two steps were repeated 2000 times to obtain the clustering results.

[0141] Based on the obtained reading clusters and the reading labels obtained by clustering, a heuristic algorithm is used to make possible adjustments to the labels of each reading. Specifically, if a reading is assigned a different label than its current one, which can obtain a better MEC score, the label will be changed. The heuristic algorithm repeatedly tries possible label adjustments for each reading. Heuristic algorithms are often prone to falling into local optimal solutions, especially when the solution space is complex or the quality of solutions varies greatly. Through its "annealing" process, simulated annealing has a higher probability of accepting inferior solutions in the early stages of the algorithm. As the "temperature" gradually decreases, the possibility of accepting inferior solutions is gradually reduced, which helps the algorithm jump out of the local optimal solution and explore a wider solution space. Through the global search capability of simulated annealing, the robustness of the heuristic algorithm to initial values ​​or random fluctuations can be improved, so that the algorithm can maintain relatively stable performance on different runs and different input data to find the best MEC score.

[0142] S52 model loss function design

[0143] The overall loss function is L = L con +β2L reg +λL AE , β2 and λ are weight parameters. The main purpose is to encourage the model to learn the characteristics of the reads more effectively so as to accurately cluster the reads from the same haplotype.

[0144] The contrastive loss is used to effectively handle the relationship between readings. The specific form is as follows:

[0145]

[0146] RetNet model output feature X 3 , the inner product matrix is ​​calculated as the similarity matrix S, where S = X 3 ·(X 3 ) T , S ij represents the similarity between the i-th read and the j-th read. The expected similarity matrix E is constructed based on the pseudo-labels obtained by cluster assignment. If reads i and j belong to the same haplotype, then E ij = 1 is considered as a positive sample pair, otherwise E ij = 0 is considered as a negative sample pair. β1 is a preset boundary value hyperparameter used to distinguish the similarity boundary between positive and negative pairs.

[0147] The loss function is constructed to encourage consistency between the similarity matrix output by the model and the expected similarity matrix, thereby indirectly minimizing the MEC score. When read i and read j are from the same haplotype, E ij =1, L con =(1-S ij ) 2 , the goal is to maximize S ij , when S ij When it is close to 1, 1-S ij The expected similarity matrix E and the correlation matrix S calculated by the model are constrained to be the same for two reads from the same haplotype. When read i and read j come from different haplotypes, if the expected similarity E ij =0, L con =max(0,s ij -β1) 2 When S ij When it is less than β1, the loss is 0. ij When it exceeds β1, the loss is (S ij -β1) 2 , encourage Si j Keep it below β1. The synchronization constraint expects the similarity matrix E and the similarity matrix S obtained by the model to behave the same for two reads from different haplotypes. The above two constraints can encourage the model to effectively learn the correlation between read pairs.

[0148] The design of the regularized loss function is based on the improvement of XHap's loss function. Since read pairs with overlapping sites can provide more predictions about the correlation of the read pairs, strong constraints should be imposed on read pairs with overlapping sites. First, a weight based on the total number of overlapping sites is assigned to each pair of reads. When there is no overlap between two read pairs, the weight is 0, avoiding the assignment of unreasonable weights to non-overlapping read pairs. Construct an m*m dynamic weight matrix W dy , specifically defined in the following form:

[0149] W dv (i, j) = overlap ij max(len i ,len j )

[0150] The overlap i，j Represents the total number of overlapping sites between read i and read j. i Represents the total number of effective sites of read i. This weight matrix is ​​used to emphasize the correlation between two read pairs, because read pairs with overlap can provide stronger evidence support for haplotype phasing.

[0151] Whether a specific sample pair is strongly positively correlated or strongly negatively correlated is reflected by calculating the correlation value of each pair of overlapping readings i and j: sim is the total number of overlapping sites where the two reads have consistent genotypes, k dissim is the total number of overlapping sites with conflicting genotypes, then the correlation between reads i and j is calculated as:

[0152]

[0153] For overlapping read pairs, if read i and read j are from the same haplotype, then C ij A positive value indicates a positive sample pair. If read i and read j come from different haplotypes, then C ij A negative value indicates a negative sample pair. If the two reads i and j have no overlapping sites, then C ij The value of C is 0. ij The calculation of the reading pair source is consistent with the similarity matrix S predicted by the model. Therefore, the regularization loss is defined as follows:

[0154]

[0155] where F ij is the value of the mask matrix at position (i, j). The mask matrix is ​​constructed based on the relationship between the overlaps between the readings. When reading i and reading j overlap, F ij The value of is 1.

[0156] The loss function of the convolutional encoder training process is described as follows:

[0157]

[0158] This loss function represents the reconstruction of the data and the original data t i where m is the number of samples in the batch.

[0159] S6 Experimental Verification

[0160] S61 Hyperparameter Settings

[0161] In the haplotype assembly experiment of simulated tetraploid potato 30× coverage short read data, Deephap was trained and optimized using the Adam optimizer to find the best parameters. The best parameter values ​​finally obtained were: learning_rate=le-4 / epoch, β1=1, β2=100, λ=0.1.

[0162] S62 evaluation index

[0163] Group the given reads into k clusters {C1, C2, ..., C k}, the corresponding MEC score can be calculated as:

[0164]

[0165] Among them, M i represents the i-th read, H j represents the jth reconstructed haplotype.

[0166] The MEC score is often used as a metric to characterize the accuracy of haplotype assembly methods. Another performance metric is the switching error rate (SWER), which is defined in diploid assembly as the proportion of positions where the phase of the reconstructed haplotype is incorrectly switched. SWER is easily generalized to polyploids to evaluate the switching error rate (VER). Among the above metrics, only MEC can be calculated without the true haplotype.

[0167] S63 Experimental Data

[0168] To tune the model and determine the hyperparameters, simulated tetraploid potato short read data with 30× coverage were used. Paired-end reads from the Illumina MiSeq technology were generated using the simulator ART. The average length of a single Illumina read was set to the maximum allowed by ART (2×250bp), and the average insert-sizes were set to 550bp with a standard deviation of 10bp. To evaluate the effect of sequencing depth on haplotypes, reads with an average coverage of 30x were obtained. A 10kbp region of potato chromosome 5 was randomly selected as the template sequence for Haplogenerator. According to the built-in log-normal model of Haplogenerator, random biallelic SNPs were introduced into each reference to generate synthetic tetraploid genomes, and the mean and scale of the logarithmic distance between SNPs were set to 3.0349 and 1.293, respectively. The reads were aligned to the reference genome using BWA-MEM. SNP sites were identified as sites where the alternative allele frequency (ie, variant allele frequency (VAF)) exceeded a predefined threshold, which was set to 0.2 in the experiment.

[0169] In order to verify the validity of this patent, this patent was tested on a real data set and compared and analyzed with five other popular haplotype assembly methods.

[0170] The Genome in a Bottle (GIAB) consortium provided the NA12878 dataset, which contains aligned PacBio SMRT whole-genome reads from a human subject with a coverage depth of 44x. We followed the benchmarking approach outlined by Wagner et al. Based on the BAM and BED files, reads covering high-confidence regions on specific chromosomes were selected, while reads covering only single SNP sites were excluded as these reads do not contain information useful for the phasing process. Specifically, 172,363 reads that aligned to chromosome 21 were selected, covering 28,719 SNPs. Finally, the SNPs matrix was constructed as input to the model.

[0171] The reference genome of potato chromosome RH89-039-165 was obtained. Paired Illumina HiSeq 2000 reads (length 2 × 100 bp) were obtained by sequencing tetraploid potato chromosome 5. Ten regions of 10 kbp in length were randomly selected from this reference genome as sample reference genomes. Reads were aligned to these reference genomes using BWA-MEM, and then SNPs were called with a VAF threshold set to 0.1.

[0172] Comparison between S64 haplotype assembly methods

[0173] This patent is compared with five other popular haplotype assembly methods, including XHap, CAECSeq, H-PoP, HapTree and HapCUT2. This patent is named DeepHap.

[0174] (1) Evaluation results of haplotype assembly using diploid real human data

[0175] Due to the large size of the read fragment matrix, we used the method in XHap to reconstruct the overlapping haplotype blocks and phased them together to obtain the complete reconstructed haplotype. Each haplotype block is 250 SNPs long, and adjacent blocks overlap by 50 SNPs. In addition, read fragments containing only one SNP site information were filtered out. The benchmark experiment results are shown in Table 1. The experimental results of HapTree cannot be obtained because the tool reports an error when running this data. In terms of continuity, DeepHap, XHap, and CAECseq can produce completely phsed haplotypes, while H-PoP and HapCUT2 produce more haplotype blocks. In terms of accuracy, CAECseq performs the worst in both MEC and SWER indicators. In terms of MEC and SWER, DeepHap's performance is basically the same as XHap, and XHap is slightly better.

[0176] Table 1. Performance comparison of DeepHap with other tools in human data

[0177]

[0178] (2) Evaluation results of haplotype assembly using tetraploid real potato data

[0179] Due to the lack of real haplotype benchmarks, the performance of each tool can only be evaluated by MEC. We benchmarked DeepHap with four other tools, XHap, CAECseq, H-PoP, and Ranbow, running the algorithm 5 times for each region's data and finally saving the haplotype with the smallest MEC. As shown in Table 2, XHap outperformed other tools in the data experiments of 10 regions.

[0180] Table 2 MEC scores of Deephap and other benchmark tools on real potato data

[0181]

[0182] The above is a detailed description of a haplotype assembly method based on deep learning provided by the present invention. For a person skilled in the art, any obvious changes made to it without departing from the essence of the present invention will constitute an infringement of the patent right of the present invention and will bear the corresponding legal liability.

Claims

1. A haplotype assembly method based on deep learning, characterized in that: The following steps are involved: S1: Data preprocessing, step S1 includes: S11: Based on the input reading data and reference genome, use BWA-MEM to align the readings with the reference genome to obtain a SAM file, extract the SNP position and genotype information based on the SAM file, and construct a SNPs matrix; S12: Since sequencing reads covering only one SNP site do not provide any information required for haplotype phasing, sequencing reads covering only one SNP site are removed from the constructed SNPs matrix; S13: One-Hot encoding of the SNPs matrix to convert the one-dimensional sequencing read allele sequence into two-dimensional read data t i ; S2: Get reading embeddings: One-Hot encoded sequencing reads are transformed into i Map to the specified lower dimension to obtain the low-dimensional potential representation h i ; S3: Extracting reading features, the step S3 comprises: S31: The learned potential representation h i Stacked together to form a reading embedding matrix, filling a virtual dimension to form a three-dimensional tensor T 0 ; S32: T 0 Input to the RetNet model for feature extraction to learn the feature representation of sequencing reads, and design constraints based on the loss function to indirectly optimize the MEC score; S4: Read count clustering: The SpectralNet model is used as the clustering layer to learn the similarity between sequencing read pairs based on the feature representation extracted by the RetNet model to complete the clustering process of sequencing reads; S5: cluster optimization and obtain result haplotypes, the step S5 comprises: S51: With the goal of optimizing the loss function, the reading clustering result is obtained by iteratively running the previous steps; S52: Use optimization algorithms to further refine the clustering results; finally, obtain haplotype consensus by majority voting based on all reads in the same cluster; each haplotype consensus result refers to one haplotype.

2. The haplotype assembly method based on deep learning according to claim 1, characterized in that: The step S2 is specifically as follows: Using the convolutional encoder from CAECSeq, the one-hot encoded sequencing reads t i Map to the specified lower dimension to obtain the low-dimensional potential representation h i ; The convolutional encoder contains three two-dimensional convolutional layers and one fully connected layer; the convolutional layer can capture local features by applying a set of learnable convolution kernels to the input data; by operating on the local area of ​​the input data, the spatial hierarchy of the data is preserved, and the output feature map can characterize the features of each part of the input, which is conducive to analyzing the two-dimensional spatial features of the sequencing read data; the stacking of three convolutional layers allows the model to learn feature hierarchies from low-level to high-level; the feature map is converted to a one-dimensional vector through the Flatten operation and input to the fully connected layer to obtain a low-dimensional potential representation h i ; The formulation of each layer of the encoder part is as follows: T 0 =t i Where "*" represents the convolution operation, and σ uses PReLU as the activation function; and Represents the weight and bias of the lth convolutional layer, W1 dense and represents the weight and bias of the last fully connected layer, h i is the 128-dimensional potential representation obtained; the sizes of the convolution kernels in each layer are (4, 5), (1, 5), (1, 3), the filter sizes are 32, 64, 128, and the step size is 1; the decoder part is formulated as follows: It is the reconstructed data obtained after the decoder; the convolutional layer design of the decoder is symmetrical with that of the encoder.

3. The haplotype assembly method based on deep learning according to claim 1, characterized in that: The step S3 further comprises S33: The RetNet model is used for feature extraction. RetNet is composed of L identical blocks stacked together, and each RetNet block contains two modules: a multi-scale retention (MSR) module and a feed-forward network (FFN) module. The core of multi-scale retention is the retention mechanism. First, the learned potential representation h i Stacked together to form a reading embedding matrix, filling a virtual dimension to form a three-dimensional tensor X 0 Input into RetNet; The L-layer model is constructed by stacking multi-scale retention modules and feed-forward network modules; this method uses X 0 As input, the L-layer RetNet model calculates the output X l : Y l =MSR(LN(X l ))+X l X l+1 =FFN(LN(Y l ))+Y l Where LN(·) is LayerNorm, FFN(X)=gelu(XW1)W2, W1 and W2 are parameter matrices; A parallel Retention mechanism is used; given an input matrix X, the parallel representation of the Retention layer is defined as follows: Mode: Q=(XW Q )⊙Θ, V=XW V Retention(X)=(QK T ☉D)V is the complex conjugate of Θ, combining causal masking and exponential decay along relative distances into a matrix D; Q n , K n is the matrix to be learned; The Retention submodule of each layer of RetNet is also divided into h heads, each head uses a different W q , W k , W v Parameters, and each head uses a different γ constant. In our model, h = 4; then for input X, the output of the MSR layer is: γ=1-2 -5-arange(0,h) head i =Retention(X,γ i ) Y=GroupNorm h (Concat(head1,...,head h )) MSR(X)=(swish(XW G )☉Y)W O Among them, swish is the activation function used to generate the gating threshold. Since each head uses a different γ parameter, the output of each head needs to be normalized and then concat.

4. The haplotype assembly method based on deep learning according to claim 1, characterized in that: The step S4 is specifically as follows: Based on the feature representation of sequencing reads obtained by the RetNet model, SpectralNet was used as the clustering module to cluster the reads into k clusters. After the optimization algorithm, the consensus sequence was calculated based on the majority vote to construct the haplotype; The entire SpectralNet training can be divided into three parts:

1. Taking features as input, a twin network is used to calculate the affinity matrix using a given distance metric; 2. Optimize the clustering objective F by enforcing orthogonality θ , maps the input data points into the eigenspace of their associated graph Laplacian matrix; 3. In the embedding space, use the k-means algorithm to learn cluster assignments; X 3 is the feature of the sequencing reads output by the last layer of RetNet model, and the twin network converts X 3 Each data point x i Mapped into a space Typically, the network is trained to minimize contrast loss. After the twin network is trained, we use it to construct an affinity matrix A as the basis for training SpectralNet. The specific calculation formula is as follows: Where c is the marginal value, usually set to 1, and σ is the standard deviation of the Gaussian distribution; SpectralNet uses a general neural network to achieve the mapping F θ , that is, y i =F θ (x i ), whose last layer has an orthogonality constraint; the model is trained by minimizing the loss LSpectralNet, which is defined as follows: Finally, for y1, ..., y m Use k-means to perform clustering and obtain clustering results.

5. The deep learning-based haplotype assembly method according to claim 1, characterized in that: The step S51 is specifically as follows: An unsupervised deep learning approach was used to extract similarities between reads and then cluster them for haplotype assembly; this was achieved through an iterative training process of a neural network. Specifically, each training epoch consisted of the following two steps: (1) In the previous step, the feature representation and cluster label of each reading can be obtained; then, the loss function is used to train the convolutional encoder model and RetNet; (2) After obtaining the optimized features, SpectralNet is used again to obtain the cluster labels of the readings; These two steps were repeated 2000 times to obtain the clustering results; The cluster assignments and read clusters obtained from model training are used as inputs to the optimization algorithm, which outputs the final optimized cluster reads. Finally, the consensus sequence is extracted as the result haplotype using the majority voting method. The overall loss function is L = L con +β2L reg +λL AE , β2 and λ are weight parameters; they are mainly used to encourage the model to learn the characteristics of the reads more effectively so as to accurately cluster the reads from the same haplotype. The contrastive loss is used to effectively handle the relationship between readings. The specific form is as follows: RetNet model output feature X 3 , the inner product matrix is ​​calculated as the similarity matrix S, where S = X 3 ·(X 3 ) T , S ij represents the similarity between the i-th reading and the j-th reading; the expected similarity matrix E is constructed based on the pseudo-labels obtained by cluster assignment. If reading i and reading j belong to the same haplotype, then E ij = 1 is considered as a positive sample pair, otherwise E ij = 0 is regarded as a negative sample pair; β1 is a preset boundary value hyperparameter used to distinguish the similarity boundary between positive and negative pairs; The loss function is constructed to encourage the consistency between the similarity matrix output by the model and the expected similarity matrix, thereby indirectly minimizing the MEC score; when reads i and j are from the same haplotype, E ij =1, L con =(1-S ij ) 2 , the goal is to maximize S ij , when S ij When it is close to 1, 1-S ij Approaching 0, thus reducing the loss; Synchronously constraining the expected similarity matrix E and the correlation matrix S calculated by the model to behave the same in two reads from the same haplotype; When read i and read j come from different haplotypes, if the expected similarity E ij =0, L con =max(0,s ij -β1) 2 ; When S ij When it is less than β1, the loss is 0. ij When it exceeds β1, the loss is (S ij -β1) 2 , encourage S ij Keep it below β1; Synchronize the expected similarity matrix E and the similarity matrix S obtained by the model to have the same performance for two reads from different haplotypes; The above two constraints can prompt the model to effectively learn the correlation between read pairs; The loss function is designed as follows: First, each pair of reads is assigned a weight based on the total number of overlapping sites. When two read pairs do not overlap, the weight is 0. Construct an m*m dynamic weight matrix W dy , specifically defined as follows: W dy (i,j)=overlap ij max(len i ,len j ) in overlap i,j indicates the total number of overlapping sites between read i and read j; len i Represents the total number of effective sites of read i. This weight matrix is ​​used to emphasize the correlation between two read pairs, because read pairs with overlap can provide stronger evidence support for haplotype phasing; Whether a specific sample pair is strongly positively correlated or strongly negatively correlated is reflected by calculating the correlation value of each pair of overlapping readings i and j: sim is the total number of overlapping sites where the two reads have consistent genotypes, k dissim is the total number of overlapping sites with conflicting genotypes, then the correlation between reads i and j is calculated as: For overlapping read pairs, if read i and read j are from the same haplotype, then C ij A positive value indicates a positive sample pair. If read i and read j come from different haplotypes, then C ij A negative value indicates a negative sample pair; if the two reads i and j have no overlapping sites, then C ij The value of C is 0; ij The calculation of the reading pair source is consistent with the similarity matrix S predicted by the model; therefore, the regularization loss is defined as follows: Among them, F ij is the value of the mask matrix at position (i, j); The loss function of the convolutional encoder training process is described as follows: This loss function represents the reconstruction of the data and the original data t i where m is the number of samples in the batch.

Images