A SNP and INDEL detection method based on deep learning and long-read sequencing

Abstract

The present invention is applicable to the field of bioinformatics, and is specifically a SNP and INDEL detection method based on deep learning and long-read sequencing, including: performing a first preprocessing on the data to be detected, generating a first superimposed image, using the first superimposed image as the input of the SNP detection model, and outputting a VCF file of the SNP; performing a second preprocessing on the data to be detected, generating a second superimposed image, using the second superimposed image as the input of the INDEL detection model, and outputting a VCF file of the INDEL; using the VcfR data processing tool based on the R language to perform variation detection and processing on the data in the obtained VCF file to obtain the variation detection results on the genome; using IGV output to display the variation detection results on the genome. The present invention uses deep learning technology and long-read sequencing technology, which can solve the low performance of traditional methods in flexibility, generalization ability, and computational efficiency, as well as the problem of undetectable areas in short-read sequencing when facing complex problems.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and specifically to a SNP and INDEL detection method based on deep learning and long-read sequencing. Background Art

[0002] Single nucleotide polymorphisms (SNPs) and insertion / deletion sequences (INDELs) are two common genetic variations in the human genome and are closely related to human genetic diseases. Using efficient SNP and INDEL detection algorithms, variations associated with specific genetic diseases can be identified in individual genomes, which helps in the diagnosis and screening of diseases.

[0003] Currently, the industry generally relies on traditional short-read sequencing methods, which are still lagging behind in flexibility, generalization ability, and computational efficiency. The computational cost required to deal with complex problems is too high, and short-read sequencing still has many difficult-to-detect areas when facing complex genomes. Summary of the invention

[0004] The purpose of an embodiment of the present invention is to provide a SNP and INDEL detection method based on deep learning and long-read sequencing. The method uses deep learning technology and long-read sequencing technology to solve the low performance of traditional related technologies in flexibility, generalization ability, and computational efficiency, as well as the technical problem of undetectable areas in short-read sequencing when facing complex problems.

[0005] To achieve the above object, the present invention provides the following technical solutions.

[0006] An embodiment of the present invention provides a SNP and INDEL detection method based on deep learning and long-read sequencing, comprising:

[0007] Obtain the data to be tested, which includes long-read sequencing data and reference genome data;

[0008] Performing a first preprocessing on the data to be detected to generate a first overlay image, using the first overlay image as an input of the SNP detection model, and outputting a VCF file of the SNP;

[0009] Performing a second preprocessing on the data to be detected to generate a second overlay image, using the second overlay image as an input of the INDEL detection model, and outputting a VCF file of the INDEL;

[0010] Using the VcfR data processing tool based on R language, the data in the obtained VCF file is subjected to variation detection and processing to obtain the variation detection results on the genome;

[0011] Use IGV output to display variant detection results on the genome.

[0012] Furthermore, the step of performing a first preprocessing on the data to be detected to generate a first superimposed image includes a process of selecting candidate sites; wherein:

[0013] Set a genomic site b with a reference base R and calculate the coverage. The calculation formula of coverage is expressed as:

[0014] ;

[0015] Where N represents the number of reads covering genomic site b; represents the long read of the i-th read segment, and L represents the length of the target region;

[0016] When the coverage of the current candidate site is lower than the coverage threshold (the threshold is set to 10), the current candidate site is removed to obtain a subset of SNP candidate sites, and the subset is used as the heterozygous SNP site set V;

[0017] Also included is a process for generating an overlay image of haplotyped SNPs, where:

[0018] For a SNP candidate site b, select sites from the site set V that share at least one read with the candidate site b and are at most 50,000 bp away from the candidate site b;

[0019] In each direction of the candidate site b, 20 sites are selected from the site set V, and Z is used to represent the set of candidate site b and potential heterozygous SNP sites near the candidate site b;

[0020] The read set covering the candidate site b is divided into four groups, where = {reads that support base B at b}, B∈{A,G,T,C}, reads that do not support any base at candidate site b are not used;

[0021] for For each read group with base B in , calculate the number of supported read groups for site t∈Z with base D∈{A,G,T,C} ;set up ,in, represents a variable, is a function that returns −1 if D is the reference basis of point t, otherwise it returns 1; thus, a 4×41×4 three-dimensional matrix M is obtained, where the entries of the three-dimensional matrix M are expressed as , the first dimension of the three-dimensional matrix M corresponds to the nucleotide type B of the B site, the second dimension corresponds to the number of t sites, and the third dimension corresponds to the nucleotide type D of the t site;

[0022] The image of the three-dimensional matrix M has read groups as rows and various base positions as columns, with 4 channels, each of which records the frequency of different bases in a given read group at a given site;

[0023] Add another channel to the image. The newly added image is a 4×41 matrix with the number of elements , if base B is the reference base for site b, then , otherwise it is 0; in this channel, there is a row of 1 as the reference basis of point b, and a row of 0 as other bases; first encode the reference bases of the base points in the set Z, and then add the encoded result as a new row to the image to obtain a superimposed image.

[0024] Furthermore, the SNP detection model is constructed based on a convolutional neural network (CNN), which includes three convolutional layers:

[0025] The first convolutional layer uses three kernels of different dimensions and combines the convolution features into a single output; one of the three kernels of different dimensions captures local information from a row, another captures local information from a column, and the last captures local information from a two-dimensional local region;

[0026] The second and third convolutional layers both use a kernel of size 2×3; the output of the third convolutional layer is flattened and used as the input to a fully connected network with dropout;

[0027] The first fully connected layer of the fully connected network is followed by two different fully connected layers of neural networks for calculating two types of probabilities; wherein, in the first neural network, the probability of each B-type nucleotide is calculated to indicate that B exists at the candidate site of the genome, and for each B-type nucleotide, there is a binary label prediction; the second neural network combines the logit output of the first neural network with the output of the first fully connected hidden layer to estimate the probability of zygosity of the candidate site.

[0028] Furthermore, the step of estimating the probability of zygosity of the candidate site comprises:

[0029] A neural network is used to calculate the probability of each nucleotide type existing at the candidate site;

[0030] If a candidate site has at least two nucleotide types with a probability greater than 0.5, it is considered heterozygous; otherwise it is considered homozygous;

[0031] For heterozygous sites, the two nucleotide types with the highest probability are selected with the heterozygous variant; for homozygous sites, only the nucleotide type with the highest probability is selected: if the nucleotide type is not the reference allele, a homozygous variant call is made; otherwise, it is homozygous-reference;

[0032] Each variant call is assigned a quality score, calculated as:

[0033] ;

[0034] In the formula, Represents the probability of the alternative allele B and is recorded as a floating point number in the QUAL field of the VCF file to indicate the chance of a false positive prediction; the larger the score, the less likely the prediction is wrong.

[0035] Furthermore, the step of outputting a VCF file of INDEL includes:

[0036] Predict SNP calls using WhatsHap;

[0037] Outputs the unphased VCF file generated by the SNP detection model and the phased VCF file generated by WhatsHap.

[0038] Furthermore, the step of performing a second preprocessing on the data to be detected and generating a second overlay image includes selecting candidate points and generating an overlay image, wherein before selecting the candidate points, a process of quality filtering the data to be detected is included, and the process is as follows:

[0039] Phred quality score is used to filter low-quality bases. The Phred quality score calculation formula is expressed as:

[0040] ;

[0041] Wherein, P is the error probability of base calling. When the quality score is lower than the threshold, the base is discarded; the threshold is set to 20;

[0042] Remove adapter sequences introduced during sequencing to avoid interference with subsequent analysis;

[0043] Remove duplicate reads generated during PCR amplification to reduce bias;

[0044] Align the reads to the reference genome and remove unaligned or low-quality alignment reads;

[0045] Remove reads that are too short or too long to ensure consistency in read length and accuracy of analysis;

[0046] Reads containing a lot of repetitive sequences or low-complexity regions were removed.

[0047] Furthermore, the step of selecting candidate points in the second preprocessing includes:

[0048] For i∈{0,1}, calculate the following three data:

[0049] Data 1: = total number of reads at stage i of site b;

[0050] Data 2: Insertion frequency = ;

[0051] Data 3: Deletion frequency = ;

[0052] B is considered an indel candidate site when the following two conditions occur:

[0053] and are greater than the specified depth threshold;

[0054] The insertion frequency is greater than the specified insertion frequency threshold, or the deletion frequency is greater than the specified deletion frequency threshold.

[0055] Furthermore, the step of generating a superimposed image in the second preprocessing includes:

[0056] For an indel candidate site b, use , and Respectively, they represent the set of all reads, reads in one phase, and reads in another phase;

[0057] set up is a reference sequence starting from site b with a length of 160 bp;

[0058] For each set , perform the following operations: For each read segment r∈S, let is a 160 bp long subsequence of the read segment starting from site b; MUSCLE was used to perform multiple sequence alignment on the following set of sequences ;set up is the re-aligned sequence, where represents the realigned reference sequence, Realignment of sequences ; Truncate all sequences of length 128 from the end; for B∈{A,G,T,C,−} and 1≤p≤128, calculate:

[0059] ;

[0060] ;

[0061] In the formula, represents a function where if If the cardinality at index p of is B, it returns 1, otherwise it returns 0; the matrix C represents the original count , p represents each symbol; M is the number of entries A 5×128 matrix;

[0062] Construct a 5×128 matrix Q whose entries are ,if At index p there is symbol B, then , otherwise 0; the matrices M and Q both have the first dimension corresponding to the symbols {A, G, T, C, −}, and the second dimension corresponding to the stacked columns of the rearranged sequence;

[0063] Construct a 5×128×2 matrix , the first channel is the matrix M−Q. The image is actually a two-dimensional matrix. What we get here is a three-dimensional matrix, which can be regarded as a stack. The second channel is the matrix Q.

[0064] Concatenate three matrices , and Together we get a 15×128×2 matrix as the input overlay image for the detection model.

[0065] Furthermore, the INDEL detection model is constructed based on a convolutional neural network (CNN), which includes three convolutional layers:

[0066] The first convolutional layer uses three kernels of different dimensions and combines the convolution features into a single output; one of the three kernels of different dimensions captures local information from a row, another captures local information from a column, and the last captures local information from a two-dimensional local area;

[0067] The second and third convolutional layers both use a kernel of size 2×3; the output of the third convolutional layer is flattened and used as the input to a fully connected network with dropout;

[0068] The first fully connected layer of the fully connected network is followed by two fully connected hidden layers to generate probability estimates for four zygosity conditions at the candidate loci: homozygous-reference, homozygous-alternative, heterozygous-reference, and heterozygous-alternative.

[0069] Furthermore, the step of using the second superimposed image as the input of the INDEL detection model and outputting a VCF file of the INDEL includes:

[0070] The probabilities of the four zygotes are calculated; if the homozygous-reference marker has the highest probability, no variant call is made; if the homozygous-alternative marker has the highest probability, a variant is called from The consensus sequence was determined from the multiple sequence alignment of and the pairwis2 global alignment algorithm with affine interval penalty of Bio-Python was used to align the candidate sites with the reference sequence;

[0071] Alternative alleles are inferred from Alignnet gaps in pairwise alignments of the two sequences;

[0072] In the case where any of the heterozygous predictions has the highest probability, use and The consensus sequence of each phase was determined separately and compared with the reference sequence. The INDEL detection results of the two phases were combined, and the final phased INDEL detection results were exported at the candidate site. The VCF file of INDEL was created based on the detection results.

[0073] Compared with the prior art, the SNP and INDEL detection method based on deep learning and long-read sequencing of the present invention has the following beneficial effects:

[0074] First, the present invention uses a more suitable preprocessing method. Specifically, the data is preprocessed using two data preprocessing methods, namely, coverage calculation and Phred quality scoring, which can better filter out low-quality sites and screen out sites with higher detection potential, and can significantly improve the accuracy of SNP and INDEL detection.

[0075] Second, compared with other long-read variant callers, the present invention can generate SNP / indel detections in complex genomic regions, and by taking advantage of long-read data in repetitive regions, it can detect SNPs / indels outside high-confidence regions that cannot be reliably detected by short-read sequencing technology, thereby providing more candidate SNP / indels sites for causal variant studies of undiagnosed diseases where short-read sequencing has not found causal candidate variants;

[0076] Third, the present invention improves the accuracy and efficiency of SNP and INDEL detection by constructing stacked images of candidate sites and using convolutional neural networks for feature extraction and classification. The present invention improves the continuity and accuracy of gene assembly by utilizing long-read sequencing data to cover a wide range of genomic regions, thereby solving the difficulties encountered when facing complex problems using ordinary short-read sequencing methods.

[0077] In summary, the present invention uses deep learning technology and long-read sequencing technology, which can effectively solve the low performance of traditional methods in flexibility, generalization ability, computational efficiency, and the problem of undetectable areas in short-read sequencing when facing complex problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present invention.

[0079] Figure 1A flow chart of a SNP and INDEL detection method based on deep learning and long-read sequencing of the present invention;

[0080] Figure 2 A schematic diagram of a process for generating a superimposed image provided by an embodiment of the present invention;

[0081] Figure 3 A schematic diagram of the flow of a convolutional neural network (CNN) provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0082] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0083] Single nucleotide polymorphisms (SNPs) and insertion / deletions (INDELs) are two common genetic variations in the human genome and are closely related to human genetic diseases. Using efficient SNP and INDEL detection algorithms, variations associated with specific genetic diseases can be identified in individual genomes, which helps diagnose and screen for diseases.

[0084] Currently, the industry generally relies on traditional short-read sequencing methods, which are still lagging behind in flexibility, generalization ability, and computational efficiency. The computational cost required to deal with complex problems is too high, and short-read sequencing still has many difficult-to-detect areas when facing complex genomes.

[0085] To solve the above problems, the present invention provides a SNP and INDEL detection method based on deep learning and long-read sequencing. The method uses the superimposed image obtained by processing long-read sequencing data and reference genome data as the input of a convolutional neural network (CNN), and uses the convolutional neural network (CNN) to generate a VCF (Variant CallFormat) file for predicting SNPs and INDELs; then the VCF file is parsed, and the variant data is processed and analyzed, and finally the IGV (Integrative Genomics Viewer) is used to view the variants on the genome.

[0086] The following is a detailed description of the specific implementation of a SNP and INDEL detection method based on deep learning and long-read sequencing of the present invention in conjunction with specific embodiments.

[0087] Please refer to Figure 1 In a preferred embodiment of the present invention, a SNP and INDEL detection method based on deep learning and long-read sequencing is provided, comprising the following steps:

[0088] S101: Acquire data to be tested, where the data to be tested includes long-read sequencing data and reference genome data;

[0089] S102: performing a first preprocessing on the data to be detected to generate a first superimposed image, using the first superimposed image as an input of the SNP detection model, and outputting a VCF file of the SNP;

[0090] S103: performing a second preprocessing on the data to be detected to generate a second superimposed image, using the second superimposed image as an input of the INDEL detection model, and outputting a VCF file of the INDEL;

[0091] S104: Using the VcfR data processing tool based on R language, perform variation detection and processing on the data in the obtained VCF file to obtain variation detection results on the genome; and use IGV output to display the variation detection results on the genome.

[0092] For further information, please refer to Figure 2 In the embodiment of the present disclosure, the step of performing a first preprocessing on the data to be detected to generate a first superimposed image includes a process of selecting candidate sites; wherein:

[0093] Set a genomic site b with a reference base R and calculate the coverage. The calculation formula for coverage is expressed as:

[0094] ;

[0095] Where N represents the number of reads covering genomic site b; represents the long read of the i-th read segment, and L represents the length of the target region;

[0096] When the coverage of the current candidate site is lower than the coverage threshold (the threshold is set to 10), the current candidate site is removed to obtain a subset of SNP candidate sites, and the subset is used as the heterozygous SNP site set V;

[0097] Furthermore, the present invention extracts long-range haplotype information from this subset of possible heterozygous SNP sites and creates an input image for all SNP candidate sites, i.e., an overlay image of haplotype SNPs for a convolutional neural network.

[0098] Among them, the subset consists of SNP candidate sites, and the allele frequency of each candidate site is calculated. The present invention specifies its coverage as 15. For heterozygous site filtering, the default range is 20 to 40, and this range can be specified by the user based on sequencing technology and read length.

[0099] For details, please refer to Figure 2 , the preferred embodiment of the present disclosure also includes a process of generating an overlay image of haplotype SNPs, wherein:

[0100] For a SNP candidate site b, a site is selected from the site set V that shares at least one read with the candidate site b and is at most 50,000 bp away from the candidate site b; it should be noted that for SNP calls in the PacBio dataset, the present invention sets this limit to 20,000 bp;

[0101] In each direction of the candidate site b, whether upstream or downstream, 20 sites are selected from the site set V. If there are less than 20 such sites, 0 can be added to the final image. Furthermore, Z is used to represent the set of potential heterozygous SNP sites near the candidate site b and the candidate site b.

[0102] The read set covering the candidate site b is divided into four groups, where = {reads that support base B at b}, B∈{A,G,T,C}, reads that do not support any base at candidate site b are not used;

[0103] for For each read group with base B in , calculate the number of supported read groups for site t∈Z with base D∈{A,G,T,C} ;set up ,in, represents a variable, is a function that returns −1 if D is the reference basis of point t, otherwise it returns 1; thus, a 4×41×4 three-dimensional matrix M is obtained, where the entries of the three-dimensional matrix M are expressed as , the first dimension of the three-dimensional matrix M corresponds to the nucleotide type B of the B site, the second dimension corresponds to the number of t sites, and the third dimension corresponds to the nucleotide type D of the t site;

[0104] In addition, the image of the three-dimensional matrix M of the present invention has read groups as rows and various base positions as columns, and has 4 channels, each of which records the frequency of different bases in a given read group at a given site;

[0105] Add another channel to the image. Here, adding a new channel to the image can be understood as superimposing an image on top of the original 4-channel image to make it a 5-channel image. The newly superimposed image is a 4×41 matrix with the number of elements being , if base B is the reference base for site b, then , otherwise it is 0; in this channel, there is a row of 1 as the reference basis of point b, and a row of 0 as other bases; first encode the base point reference bases in the set Z, and then add the encoded result to the image as a new row to obtain a superimposed image; wherein, the new row requires first encoding the base point reference bases in Z, and then adding the encoded result to the image as a new row.

[0106] Please refer to Figure 3 In the disclosed embodiment, the SNP detection model is constructed based on a convolutional neural network (CNN), which includes three convolutional layers: the first convolutional layer uses three kernels of different dimensions and combines the convolutional features into a single output; one of the three kernels of different dimensions captures local information from a row, another captures local information from a column, and the last captures local information from a two-dimensional local area;

[0107] The second and third convolutional layers both use a kernel of size 2×3; the output of the third convolutional layer is flattened and used as the input to a fully connected network with dropout;

[0108] The first fully connected layer of the fully connected network is followed by two different fully connected layers of neural networks for calculating two types of probabilities; wherein, in the first neural network, the probability of each B-type nucleotide is calculated to indicate that B exists at the candidate site of the genome, and for each B-type nucleotide, there is a binary label prediction; the second neural network combines the logit output of the first neural network with the output of the first fully connected hidden layer to estimate the probability of zygosity of the candidate site.

[0109] Furthermore, in an embodiment of the present invention, the step of estimating the probability of zygosity of the candidate site includes:

[0110] A neural network is used to calculate the probability of each nucleotide type existing at the candidate site;

[0111] If a candidate site has at least two nucleotide types with a probability greater than 0.5, it is considered heterozygous; otherwise, it is considered homozygous;

[0112] For heterozygous sites, the two nucleotide types with the highest probability are selected with the heterozygous variant; for homozygous sites, only the nucleotide type with the highest probability is selected: if the nucleotide type is not the reference allele, a homozygous variant call is made; otherwise, it is homozygous-reference;

[0113] Each variant call is assigned a quality score, calculated as:

[0114] ;

[0115] In the formula, Represents the probability of the alternative allele B and is recorded as a floating point number in the QUAL field of the VCF file to indicate the chance of a false positive prediction; the larger the score, the less likely the prediction is wrong.

[0116] Furthermore, the step of outputting a VCF file of INDEL includes:

[0117] Predict SNP calls using WhatsHap;

[0118] Outputs the unphased VCF file generated by the SNP detection model and the phased VCF file generated by WhatsHap.

[0119] Please continue to refer to Figure 2 In a preferred embodiment of the present disclosure, the step of performing a second preprocessing on the data to be detected and generating a second overlay image includes selecting candidate points and generating an overlay image, wherein before selecting the candidate points, a process of quality filtering the data to be detected is included, and the process is as follows:

[0120] Phred quality score is used to filter low-quality bases. The Phred quality score calculation formula is expressed as:

[0121] ;

[0122] Wherein, P is the error probability of base calling. When the quality score is lower than the threshold, the base is discarded; the threshold is set to 20;

[0123] Remove the adapter sequences introduced during sequencing to avoid interference with subsequent analysis; remove duplicate reads generated during PCR amplification to reduce bias; align the reads with the reference genome and remove reads that are not aligned or have low alignment quality; remove reads that are too short or too long to ensure consistency in read length and accuracy of analysis; remove reads containing a large number of repetitive sequences or low-complexity regions; and then select candidate sites. The step of selecting candidate sites in the second preprocessing includes:

[0124] For i∈{0,1}, calculate the following three data:

[0125] Data 1: = total number of reads at stage i of site b;

[0126] Data 2: Insertion frequency = ;

[0127] Data 3: Deletion frequency = ;

[0128] B is considered an indel candidate site when the following two conditions occur:

[0129] and are greater than the specified depth threshold;

[0130] The insertion frequency is greater than the specified insertion frequency threshold, or the deletion frequency is greater than the specified deletion frequency threshold.

[0131] Furthermore, the step of generating a superimposed image in the second preprocessing includes:

[0132] For an indel candidate site b, use , and Respectively, they represent the set of all reads, reads in one phase, and reads in another phase;

[0133] set up is a reference sequence starting from site b with a length of 160 bp;

[0134] For each set , perform the following operations: For each read segment r∈S, let is a 160 bp long subsequence of the read segment starting from site b; MUSCLE was used to perform multiple sequence alignment on the following set of sequences ;set up is the re-aligned sequence, where represents the realigned reference sequence, Realignment of sequences ; Truncate all sequences of length 128 from the end; for B∈{A,G,T,C,−} and 1≤p≤128, calculate:

[0135] ;

[0136] ;

[0137] In the formula, represents a function where if If the cardinality at index p of is B, it returns 1, otherwise it returns 0; the matrix C represents the original count , p represents each symbol; M is the number of entries A 5×128 matrix;

[0138] Furthermore, a 5×128 matrix Q is constructed, whose entries are ,if At index p there is symbol B, then , otherwise 0; the matrices M and Q both have the first dimension corresponding to the symbols {A, G, T, C, −}, and the second dimension corresponding to the stacked columns of the rearranged sequence;

[0139] Further, construct a 5×128×2 matrix , the first channel is the matrix M−Q. The image is actually a two-dimensional matrix. What we get here is a three-dimensional matrix, which can be regarded as a stack. The second channel is the matrix Q.

[0140] Concatenate three matrices , and Together we get a 15×128×2 matrix as the input overlay image for the detection model.

[0141] Therefore, the present invention constructs a stacked image of candidate sites and uses a convolutional neural network to extract and classify features to improve the accuracy and efficiency of SNP and INDEL detection; the present invention uses long-read sequencing data to cover a wide range of genomic regions to improve the continuity and accuracy of gene assembly, thereby solving the difficulties encountered when facing complex problems using ordinary short-read sequencing methods.

[0142] For further information, please refer to Figure 3 In the disclosed embodiment, the INDEL detection model is constructed based on a convolutional neural network (CNN), including three convolutional layers:

[0143] The first convolutional layer uses three kernels of different dimensions and combines the convolution features into a single output; one of the three kernels of different dimensions captures local information from a row, another captures local information from a column, and the last captures local information from a two-dimensional local region;

[0144] The second and third convolutional layers both use a kernel of size 2×3; the output of the third convolutional layer is flattened and used as the input to a fully connected network with dropout;

[0145] The above is similar to the structure of the SNP detection model. Furthermore, the difference between the SNP detection model and the INDEL detection model lies in the following parts: the first fully connected layer of the fully connected network of the INDEL detection model is followed by two fully connected hidden layers, which are used to generate probability estimates of four zygotic conditions at the candidate site: homozygous-reference, homozygous-alternative, heterozygous-reference and heterozygous-alternative.

[0146] Further, the step of using the second superimposed image as the input of the INDEL detection model and outputting the VCF file of the INDEL includes: calculating the probabilities of the four zygotes; if the homozygous-reference marker has the highest probability, no variant call is performed; if the homozygous-alternative label has the highest probability, The consensus sequence was determined from a multiple sequence alignment and aligned to the reference sequence at candidate sites using the Bio-Python pairwise2 global alignment algorithm with an affine gap penalty; alternative alleles were inferred from Alignmnet gaps in the pairwise alignment of the two sequences.

[0147] In the case where any of the heterozygous predictions has the highest probability, use and The consensus sequence of each phase was determined separately and compared with the reference sequence. The INDEL detection results of the two phases were combined, and the final phased INDEL detection results were exported at the candidate site. The VCF file of INDEL was created based on the detection results.

[0148] Compared with other long-read variant callers, the present invention can generate SNP / indel detections in complex genomic regions, and by taking advantage of long-read data in repetitive regions, it is able to detect SNPs / indels outside high-confidence regions that cannot be reliably detected by short-read sequencing technology, thereby providing more candidate SNP / indels sites for causal variant studies of undiagnosed diseases where causal candidate variants are not found by short-read sequencing.

[0149] The present invention uses two authoritative public data sets, the Nanopore long-read data set and the PacBio long-read data set, to evaluate the performance of the present invention.

[0150] The present invention is more competitive than other long-read variant callers and can generate SNP / indel detections in complex genomic regions; taking advantage of long-read data in repetitive regions, the present invention is able to detect SNP / indels outside high-confidence regions that cannot be reliably detected by short-read sequencing technology, thereby providing more candidate SNP / indels sites for causal variation studies of undiagnosed diseases where no causal candidate variation was found by short-read sequencing.

[0151] The present invention uses rescaled statistics to generate stacking candidate sites, and the rescaled statistics are independent of the coverage of the test genome, so it can handle a test data set with a coverage different from that of the training data set and has strong generalization performance.

[0152] Overall, the present invention enables the detection of genetic variation in genomic regions that were previously inaccessible to genome sequencing and may facilitate the use of long-read sequencing for disease variant discovery in human genetic studies;

[0153] It should be noted that the device embodiments described above are merely schematic, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. In addition, in the drawings of the system embodiments provided by the present invention, the connection relationship between the modules indicates that there is a communication connection between them, which may be specifically implemented as one or more communication buses or signal lines. A person of ordinary skill in the art may understand and implement it without paying any creative effort.

[0154] Technicians in the field can clearly understand that for the sake of convenience and conciseness, the specific working process of the system described above can refer to the corresponding process in the aforementioned SNP and INDEL detection method embodiment based on deep learning and long-read sequencing, and will not be repeated here.

[0155] Another preferred embodiment of the present invention provides a computer device, comprising a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor, wherein when the processor executes the computer program, a SNP and INDEL detection method based on deep learning and long-read sequencing as described in any one of the above embodiments is implemented.

[0156] The computer device may be a computing device such as a desktop computer, a notebook, a PDA, a cloud server, etc.; the identification device may include but is not limited to a processor and a memory.

[0157] The processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor, etc. The processor is the control center of the terminal device, and uses various interfaces and lines to connect various parts of the entire terminal device.

[0158] The memory can be used to store the computer program, and the processor realizes various functions of the terminal device by running or executing the computer program stored in the memory and calling the data stored in the memory. The memory can mainly include a program storage area and a data storage area, wherein the program storage area can store an operating system, at least one application required for a function, etc.; the data storage area can store data created according to the use of the mobile phone, etc. In addition, the memory can include a high-speed random access memory, and can also include a non-volatile memory, such as a hard disk, a memory, a plug-in hard disk, a smart memory card (Smart Media Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), at least one disk storage device, a flash memory device, or other volatile solid-state storage devices.

[0159] Another preferred embodiment of the present invention provides a storage medium, which is a computer-readable storage medium, and the computer program is stored in the computer-readable storage medium. When the computer program is executed by a processor, it can implement the steps of each of the above-mentioned SNP and INDEL detection method embodiments based on deep learning and long-read sequencing.

[0160] The computer program includes computer program code, which may be in source code form, object code form, executable file or some intermediate form, etc. The computer readable medium may include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard disk, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electric carrier signal, telecommunication signal and software distribution medium, etc.

[0161] The above is a preferred embodiment of the present invention. It should be pointed out that a person skilled in the art can make several improvements and modifications without departing from the principle of the present invention. These improvements and modifications are also considered to be within the scope of protection of the present invention.

Claims

1. A SNP and INDEL detection method based on deep learning and long-read sequencing, characterized in that: include: Obtain the data to be tested, which includes long-read sequencing data and reference genome data; Performing a first preprocessing on the data to be detected to generate a first overlay image, using the first overlay image as an input of the SNP detection model, and outputting a VCF file of the SNP; Performing a second preprocessing on the data to be detected to generate a second overlay image, using the second overlay image as an input of the INDEL detection model, and outputting a VCF file of the INDEL; Using the VcfR data processing tool based on R language, the data in the obtained VCF file is subjected to variation detection and processing to obtain the variation detection results on the genome; Use IGV output to display variant detection results on the genome; The step of performing a first preprocessing on the data to be detected to generate a first superimposed image includes a process of selecting candidate sites; wherein: Set a genomic site b with a reference base R and calculate the coverage. The calculation formula of coverage is expressed as: ; Where N represents the number of reads covering genomic site b; represents the long read of the i-th read segment, and L represents the length of the target region; When the coverage of the current candidate site is lower than the coverage threshold, the current candidate site is removed to obtain a subset of SNP candidate sites, and the subset is used as the heterozygous SNP site set V; Also included is a process for generating an overlay image of haplotyped SNPs, where: For a SNP candidate site b, select sites from the site set V that share at least one read with the candidate site b and are at most 50,000 bp away from the candidate site b; In each direction of the candidate site b, 20 sites are selected from the site set V, and Z is used to represent the set of candidate site b and potential heterozygous SNP sites near the candidate site b; The read set covering the candidate site b is divided into four groups, where = {reads that support base B at b}, B∈{A,G,T,C}, reads that do not support any base at candidate site b are not used; for For each read group with base B in , calculate the number of supported read groups for site t∈Z with base D∈{A,G,T,C} ;set up ,in, represents a variable, is a function that returns −1 if D is the reference basis of point t, otherwise it returns 1; thus, a 4×41×4 three-dimensional matrix M is obtained, where the entries of the three-dimensional matrix M are expressed as , the first dimension of the three-dimensional matrix M corresponds to the nucleotide type B of the B site, the second dimension corresponds to the number of t sites, and the third dimension corresponds to the nucleotide type D of the t site; The image of the three-dimensional matrix M has read groups as rows and various base positions as columns, with 4 channels, each of which records the frequency of different bases in a given read group at a given site; Add another channel to the image. The newly added image is a 4×41 matrix with the number of elements , if base B is the reference base for site b, then , otherwise it is 0; in this channel, there is a row of 1 as the reference basis of point b, and a row of 0 as other bases; first encode the reference bases of the base points in the set Z, and then add the encoded result as a new row to the image to obtain a superimposed image.

2. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 1, characterized in that: The SNP detection model is built based on a convolutional neural network (CNN), which includes three convolutional layers: The first convolutional layer uses three kernels of different dimensions and combines the convolution features into a single output; one of the three kernels of different dimensions captures local information from a row, another captures local information from a column, and the last captures local information from a two-dimensional local area; The second and third convolutional layers both use a kernel of size 2×3; the output of the third convolutional layer is flattened and used as the input to a fully connected network with dropout; The first fully connected layer of the fully connected network is followed by two different fully connected layers of neural networks for calculating two types of probabilities; wherein, in the first neural network, the probability of each B-type nucleotide is calculated to indicate that B exists at the candidate site of the genome, and for each B-type nucleotide, there is a binary label prediction; the second neural network combines the logit output of the first neural network with the output of the first fully connected hidden layer to estimate the probability of zygosity of the candidate site.

3. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 2, characterized in that: The step of estimating the probability of zygosity of the candidate site comprises: A neural network is used to calculate the probability of each nucleotide type existing at the candidate site; If a candidate site has at least two nucleotide types with a probability greater than 0.5, it is considered heterozygous; otherwise, it is considered homozygous; For heterozygous sites, the two nucleotide types with the highest probability were selected with the heterozygous variant; for homozygous sites, only the nucleotide type with the highest probability was selected; if the nucleotide type was not the reference allele, a homozygous variant call was made; otherwise, it was homozygous-reference; Each variant call is assigned a quality score, calculated as: ； In the formula, Represents the probability of the alternative allele B and is recorded as a floating point number in the QUAL field of the VCF file to indicate the chance of a false positive prediction; the larger the score, the less likely the prediction is wrong.

4. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 3, characterized in that: The step of outputting a VCF file of INDEL comprises: Predict SNP calls using WhatsHap; Outputs the unphased VCF file generated by the SNP detection model and the phased VCF file generated by WhatsHap.

5. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 4, characterized in that: The step of performing a second preprocessing on the data to be detected and generating a second overlay image includes selecting candidate points and generating an overlay image, wherein before selecting the candidate points, a process of performing quality filtering on the data to be detected includes: Phred quality score is used to filter low-quality bases. The Phred quality score calculation formula is expressed as: ； Wherein, P is the error probability of base calling. When the quality score is lower than the threshold, the base is discarded; the threshold is set to 20; Remove adapter sequences introduced during sequencing; remove duplicate reads generated during PCR amplification; align reads to the reference genome and remove reads that are not aligned or have low alignment quality; remove reads that are too short or too long; remove reads that contain a large number of repetitive sequences or low-complexity regions.

6. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 5, characterized in that: The step of selecting candidate points in the second preprocessing includes: For i∈{0,1}, calculate the following three data: Data 1: = total number of reads at stage i of site b; Data 2: Insertion frequency = ; Data 3: Deletion frequency = ; B is considered an indel candidate site when the following two conditions occur: and are greater than the specified depth threshold; The insertion frequency is greater than the specified insertion frequency threshold, or the deletion frequency is greater than the specified deletion frequency threshold.

7. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 6, characterized in that: The step of generating a superimposed image in the second preprocessing includes: For an indel candidate site b, use , and Respectively, they represent the set of all reads, reads in one phase, and reads in another phase; set up is a reference sequence starting from site b with a length of 160 bp; For each set , perform the following operations: For each read segment r∈S, let is a 160 bp long subsequence of the read segment starting from site b; MUSCLE was used to perform multiple sequence alignment on the following set of sequences ; set up is the re-aligned sequence, where represents the realigned reference sequence, Realignment of sequences ; Truncate all sequences of length 128 from the end; For B∈{A,G,T,C,−} and 1≤p≤128, calculate: ; ; In the formula, represents a function where if If the cardinality at index p of is B, it returns 1, otherwise it returns 0; the matrix C represents the original count , p represents each symbol; M is the number of entries A 5×128 matrix; Construct a 5×128 matrix Q whose entries are ,if At index p there is symbol B, then , otherwise 0; the matrices M and Q both have the first dimension corresponding to the symbols {A, G, T, C, −}, and the second dimension corresponding to the stacked columns of the rearranged sequence; Construct a 5×128×2 matrix , the first channel is the matrix M−Q. The image is actually a two-dimensional matrix. What we get here is a three-dimensional matrix, which can be regarded as a stack. The second channel is the matrix Q. Concatenate three matrices , and Together we get a 15×128×2 matrix as the input overlay image for the detection model.

8. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 7, characterized in that: The INDEL detection model is built based on a convolutional neural network (CNN), which includes three convolutional layers: The first convolutional layer uses three kernels of different dimensions and combines the convolution features into a single output; one of the three kernels of different dimensions captures local information from a row, another captures local information from a column, and the last captures local information from a two-dimensional local region; The second and third convolutional layers both use a kernel of size 2×3; the output of the third convolutional layer is flattened and used as the input to a fully connected network with dropout; The first fully connected layer of the fully connected network is followed by two fully connected hidden layers to generate probability estimates for four zygosity conditions at the candidate site: homozygous-reference, homozygous-alternative, heterozygous-reference, and heterozygous-alternative.

9. The SNP and INDEL detection method based on deep learning and long-read sequencing according to claim 8, characterized in that: The step of using the second superimposed image as the input of the INDEL detection model and outputting a VCF file of the INDEL comprises: The probabilities of the four zygotes are calculated; if the homozygous-reference marker has the highest probability, no variant call is made; if the homozygous-alternative marker has the highest probability, a variant is called from The consensus sequence was determined from the multiple sequence alignment of and the pairwis2 global alignment algorithm with affine interval penalty of Bio-Python was used to align the candidate sites with the reference sequence; Alternative alleles are inferred from Alignnet gaps in pairwise alignments of the two sequences; In the case where any of the heterozygous predictions has the highest probability, use and The consensus sequence of each phase was determined separately and compared with the reference sequence. The INDEL detection results of the two phases were combined, and the final phased INDEL detection results were exported at the candidate site. The VCF file of INDEL was created based on the detection results.

Images