Method for filtering sequencing data of living body

Abstract

The invention belongs to the technical field of bioinformatics, and particularly relates to a method for filtering sequencing data of organisms, which comprises the following steps of: clustering the sequencing data of the organisms according to the content of a characteristic fragment set and species characteristics to obtain a plurality of read length groups; marking a target read length group with a high-credibility read length sequence from the plurality of read length groups, wherein the high-credibility read length sequence is a sequence of which the content of the feature fragment set is greater than a threshold value; and repeating the cycle, and screening out the read length sequence with the occurrence frequency higher than the set value in all the target read length groups. The method does not depend on a large reference genome database, genome effective extraction is performed on offspring sequencing data by using family genetic characteristics, pollution sequences are removed, and a symbiotic genome is separated; and multiple times of loop iteration clustering are adopted, so that the accuracy of sequence classification can be improved as much as possible, and the purpose of efficiently removing polluted sequences is achieved.

Description

Technical Field

[0001] The invention belongs to the technical field of bioinformatics, and in particular relates to a method for filtering sequencing data of an organism. Background Art

[0002] At present, contamination is widely distributed in gene sequencing data, including contaminated sequences caused by various reasons such as random sequencing errors, artificial mixed contamination, and symbiotic parasites, which seriously weaken the credibility of downstream bioinformatics analysis. Therefore, the development of efficient methods for separation, identification and filtering of sequencing data is of great significance to restore the real life information conveyed by sequencing data and clarify fine genetic variations and biological interactions. However, there are limited technical methods for filtering sequencing data contamination at home and abroad, and most of them are concentrated on sequence similarity comparison and broad-spectrum sterilization technology based on second-generation high-throughput sequencing data. These traditional methods often rely on comparison with existing public nucleic acid databases to remove "contaminated" sequences. Although they can remove a certain amount of contaminated sequences, they are limited by sequencing technology, known species information and genome integrity, and the efficiency of contaminated data filtering needs to be improved.

[0003] In recent years, with the continuous development of single-molecule sequencing technology, the sequencing data it produces has the advantage of long read length, which is expected to accelerate the study of the genomics of new species and improve the quality of the genome sequences of existing species. However, single-molecule sequencing technology also faces defects such as high sequencing error rates. Traditional decontamination methods rely on reference genome sequences or large databases for accurate sequence alignment, so they cannot be directly applied to single-molecule sequencing data with high error rates, which to some extent hinders the research and application of single-molecule sequencing technology in genomics. To address this problem, the early algorithm model was based on algorithm optimization such as interval sampling to weaken the chain effect of single base errors and tolerate high read length error rates, but the improvement in feature recognition accuracy was limited. Summary of the invention

[0004] The present invention provides a method for filtering sequencing data of an organism, which realizes rapid filtering without relying on a public nucleic acid database and can accurately identify the sequence of the organism.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is:

[0006] A first aspect of the present invention provides a method for filtering sequencing data of an organism, comprising the following steps:

[0007] S110: clustering the sequencing data of the organism according to the content and species characteristics of the characteristic fragment set to obtain a plurality of read length groups, wherein the characteristic fragment set includes at least one of a parental common characteristic fragment set, a paternal unique characteristic fragment set, and a maternal unique characteristic fragment set;

[0008] S120: marking a target read length group having a high-confidence read length sequence from the plurality of read length groups, wherein the high-confidence read length sequence is a sequence whose content of the characteristic fragment set is above a threshold value;

[0009] S130: Repeat the cycle of S110 and S120 to screen out read length sequences whose occurrence frequency in all target read length groups is higher than a set value.

[0010] In some embodiments of the present invention, the characteristic fragment set is obtained by performing set operations based on the paternal Kmer sequence set and the maternal Kmer sequence set.

[0011] In some embodiments of the present invention, the characteristic fragment set obtained by performing set operations based on the paternal Kmer sequence set and the maternal Kmer sequence set includes:

[0012] Based on the sequencing data of the paternal and maternal copies of the organism, a paternal Kmer sequence set and a maternal Kmer sequence set are obtained;

[0013] A set operation is performed based on the paternal Kmer sequence set and the maternal Kmer sequence set to obtain a characteristic fragment set.

[0014] In some embodiments of the present invention, the set operation includes taking intersection and difference.

[0015] In some embodiments of the present invention, the characteristic fragment set is a Kmer sequence set.

[0016] In some embodiments of the present invention, the length of the Kmer sequence in the Kmer sequence set is 10 to 100 bases.

[0017] In some embodiments of the present invention, before the set operation, the paternal Kmer sequence set and the maternal Kmer sequence set are converted into a paternal Strobemer sequence set and a maternal Strobemer sequence set using the Strobemer method; the characteristic fragment set is a Strobemer sequence set.

[0018] In some embodiments of the present invention, the strobemer sequence in the strobemer sequence set is based on a sliding window of variable size, each window extracts an l-mer, and connects n consecutive l-mers head to tail to obtain n*l-mer; wherein l is 10 to 100, n is 2 to 10, and the variable size has a minimum value w min and the maximum value w max , w min 5~20,w max is 10 to 50, and w min <wmax .

[0019] In some embodiments of the invention, the Strobemer method comprises at least one of minstrobes, hybridstrobes, randstrobes, mixedstrobes, altstrobes, and multistrobes.

[0020] In some embodiments of the present invention, the sequencing data of the father and mother of the organism are second-generation sequencing data.

[0021] In some embodiments of the present invention, the sequencing data of the organism is third-generation sequencing data.

[0022] In some embodiments of the present invention, the third-generation sequencing data is any one of PacBio HiFi, PacBio CLR, ONT, and Cyclone WT sequencing data.

[0023] In some embodiments of the present invention, in S120, the sequence whose content of the characteristic fragment set is above the threshold is a sequence whose content of the father's specific characteristic fragment set is above the first threshold, whose content of the mother's specific characteristic fragment set is above the second threshold, and whose content of the parent-common characteristic fragment set is above the third threshold.

[0024] In some embodiments of the present invention, the first threshold, the second threshold and the third threshold are independently 0.01-1%.

[0025] In some embodiments of the present invention, the first threshold is independently 0.05-0.5%, the second threshold is independently 0.05-0.5%, and the third threshold is independently 0.2-0.8%.

[0026] In some embodiments of the present invention, the first threshold is independently 0.05-0.2%, the second threshold is independently 0.05-0.2%, and the third threshold is independently 0.4-0.6%.

[0027] In some embodiments of the present invention, the species characteristics include at least one of GC content and frequency of oligonucleotide markers.

[0028] In some embodiments of the present invention, the length of the oligonucleotide tag is 2 to 10.

[0029] In some embodiments of the present invention, the length of the oligonucleotide tag is 2 to 4.

[0030] In some embodiments of the invention, the length of the oligonucleotide tag is 3.

[0031] In some embodiments of the present invention, among the frequencies of the oligonucleotide markers, the frequencies of the complementary oligonucleotides are combined and calculated.

[0032] In some embodiments of the present invention, in S130, the number of cycles is 5 to 20 times.

[0033] In some embodiments of the present invention, the set value is 60-99%.

[0034] A second aspect of the present invention provides a method for assembling sequencing data of an organism, comprising the following steps:

[0035] Filter according to the aforementioned filtering method, and then assemble the filtered read sequences.

[0036] In some embodiments of the present invention, the assembled tool is selected from Canu and Flye.

[0037] In some embodiments of the present invention, the method further includes classifying the read sequences into paternal-specific read sequences and maternal-specific read sequences according to the paternal-specific characteristic fragment set and the maternal-specific characteristic fragment set, and assembling them separately.

[0038] According to a third aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores computer-executable instructions, and the computer-executable instructions are used to enable a computer to execute the aforementioned filtering method or assembling method.

[0039] A fourth aspect of the present invention provides an electronic device, comprising a processor and a memory, wherein the memory stores a computer program executable on the processor, and the processor implements the aforementioned filtering method or assembly method when running the computer program.

[0040] A fifth aspect of the present application provides a sequencing data filtering system, the sequencing data filtering system comprising:

[0041] A clustering module, wherein the clustering module is used to cluster the sequencing data of the organism according to the content and species characteristics of the characteristic fragment set to obtain a plurality of read length groups, wherein the characteristic fragment set includes at least one of a parental common characteristic fragment set, a paternal unique characteristic fragment set, and a maternal unique characteristic fragment set;

[0042] A marking module, the marking module is used to mark a target read length group having a high-confidence read length sequence from the plurality of read length groups, the high-confidence read length sequence being a sequence whose content of the characteristic fragment set is above a threshold value;

[0043] A screening module is used to screen out read sequences whose occurrence frequency in all target read length groups is higher than a set value after the clustering module and the marking module are repeatedly cycled.

[0044] A sixth aspect of the present application provides a sequencing data assembly system, comprising:

[0045] The aforementioned filtration system;

[0046] An assembly module is used to assemble the read length sequences screened by the filtering system.

[0047] The beneficial effects of the present invention are:

[0048] (1) The present invention does not rely on a large reference genome database, but uses pedigree genetic characteristics to effectively extract genomes from progeny sequencing data, remove contaminating sequences, and separate symbiotic genomes.

[0049] (2) The present invention uses fragment units of Kmer sequences or error-tolerant fragment units of Strobemer sequences to effectively obtain characteristic information of single-molecule data with a high error rate.

[0050] (3) The present invention uses multiple cycles of iterative clustering to maximize the accuracy of sequence classification, thereby achieving the purpose of efficiently removing contaminated sequences.

[0051] (4) The present invention uses filtered target species data for higher-precision assembly, and can also be applied to symbiont samples for metagenomic assembly, which can effectively restore the original appearance of the host and symbiont genomes. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 It is a flowchart of the method for filtering sequencing data in some embodiments of the present application.

[0053] Figure 2 It is the statistical distribution of GC content of long read data of three platforms simulated in some embodiments of the present application.

[0054] Figure 3 It is the statistical distribution of GC content of long read data of three platforms actually sequenced in some embodiments of the present application. DETAILED DESCRIPTION

[0055] In the description of the present invention, the terms "first", "second", and "third" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Thus, a feature defined as "first", "second", and "third" may explicitly or implicitly include at least one of the features. In the description of the present invention, the meaning of "plurality" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

[0056] A first aspect of the present invention provides a method for filtering sequencing data of an organism.

[0057] Among them, an organism refers to an organism developed from a fertilized egg, which is usually produced by hybridization of a sexually reproduced father and mother. Relative to the father and mother, it is called an offspring. The offspring will inherit half of the chromosomes of the father and mother, so the genome of the offspring contains a genome sequence with genetic information of a part of the father and mother. In some embodiments, the organism includes sexually reproduced protists and fungi, as well as animals and plants. In some embodiments, sexually reproduced plants include higher plants, such as mosses, ferns, gymnosperms and angiosperms. In some embodiments, sexually reproduced animals include chordates, such as vertebrates, including at least one of fish, amphibians, reptiles, birds, mammals, etc. In some embodiments, mammals include but are not limited to Carnivora, Perissodactyla, Artiodactyla, Rodentia, Lagomorpha, Primates, etc., such as dogs, cats, pigs, sheep, cattle, horses, rabbits, mice, monkeys, orangutans, gorillas, chimpanzees, humans, etc. At least one of them.

[0058] Wherein, the organism extracts genetic material from the sample and performs sequencing to obtain corresponding sequencing data, and the sample type includes at least one of environmental samples and biological samples. In some embodiments, environmental samples include but are not limited to at least one of water bodies (such as domestic water, industrial water, medical water, agricultural water and other types of sewage and wastewater, or rivers, ocean water bodies, etc.), air, soil, compost, sludge (such as river sludge, wastewater sedimentation tank sludge, etc.), volcanic ash, frozen soil, and food (such as solid food, fluid food, beverages, etc.). In some embodiments, biological samples include but are not limited to at least one of body fluids (such as blood, tissue fluid, lymph, cerebrospinal fluid, urine, sweat, sputum, saliva, gastric juice, intestinal juice, pancreatic juice, bile, prostatic fluid, vaginal secretions, semen, serous cavity effusion, joint cavity effusion, bronchoalveolar lavage fluid, amniotic fluid, etc.), skin, feces, intestinal contents, and histological samples (such as samples obtained by surgery, endoscopy or percutaneous puncture biopsy). In some embodiments, the genetic material is at least one of DNA and RNA. Since mitochondrial DNA or chloroplast DNA is usually inherited from the mother, in the embodiments of the present application, the DNA is usually chromosomal DNA.

[0059] In the process of extracting, separating, and sequencing samples of the organism, due to factors such as the mixing of symbiotic or parasitic organisms, accidental mixing of exogenous cells in the laboratory, unclean sequencing vessels, and sequencing errors, the sequencing data obtained by sequencing the genetic material extracted from the sample may contain, in addition to the sequences belonging to the organism itself, stable or random contamination sequences caused by one or several of the above factors, such as symbiotic organism sequences, parasitic organism sequences, etc.

[0060] Sequencing data can be obtained by any of the methods of first-generation sequencing, second-generation sequencing, and third-generation sequencing. Among them, first-generation sequencing includes Maxam-Gilbert sequencing technology, Sanger dideoxy sequencing technology, fluorescent automatic sequencing technology, and hybrid sequencing technology, etc., second-generation sequencing includes 454Roche GS FLX, Illumina, SOLiD, Ion Torrent, BGISEQ, etc., and third-generation sequencing includes HeliScope, PacBio HiFi, PacBio CLR, ONT, Cyclone WT, etc. Third-generation sequencing is also called single-molecule sequencing, which can read nucleotide sequences at the single-molecule level, and has the advantages of longer read length and faster sequencing speed. Therefore, in some embodiments of the present application, the genetic material extracted from the organism is sequenced by single-molecule sequencing to obtain corresponding sequencing data.

[0061] refer to Figure 1 , the method for filtering sequencing data of an organism includes the following steps S110 to S130.

[0062] S110: Clustering the sequencing data of the organism according to the content of the characteristic fragment set and the species characteristics to obtain a plurality of read length groups, wherein the characteristic fragment set includes at least one of a parental common characteristic fragment set, a paternal unique characteristic fragment set and a maternal unique characteristic fragment set.

[0063] The organism in the embodiment of the present application belongs to the family genetic model, and the offspring organism has a genome sequence with genetic information of each part of the father and mother. Because the father and mother are from the same species, there are more homologous sequences in the father's genome and the mother's genome; but due to gender differences and genetic diversity, there are also certain heterogeneous sequence sites between the father and the mother. In the genetic process, the offspring will obtain half of the chromosomes of the father and the mother through inheritance, and these chromosomes carry the characteristics of the homologous sequences, heterogeneous sequences and germline mutations of the father and the mother to varying degrees, which can be specifically embodied by the specific heterozygous sites of the father and the mother's chromosomes. In the process of extracting, separating, sequencing and other processes of the sample of the organism, due to the mixing of symbiotic or parasitic organisms, the mixing of foreign cells in the laboratory, the uncleanness of the sequencing vessel and other factors, the sequencing sequence of the organism measured may contain stable or random contaminants. These contaminants can be divided into the following types according to their different coexistence relationships with the three individuals in the family of the organism: those that only exist in the offspring sequencing batch (accidental sampling or experimental contamination), those that only exist in the offspring and the paternal individual (accidental contamination), those that only exist in the offspring and the maternal individual (accidental contamination), and those that are shared by all individuals (relatively stable symbiotic and parasitic species). Therefore, by performing set operations on the paternal, maternal, and offspring data containing different contamination sequences, the offspring organism genome sequencing sequences with strong sequence characteristics can be extracted. In theory, the paternal and maternal specific sequence sets can mark the homologous chromosome heterozygous intervals inherited by the offspring, eliminating contamination that only exists in the offspring and is shared; while the paternal and maternal common sequence sets can identify the species conservative intervals in the offspring individuals, eliminating accidental contamination that only exists in the father and son or the mother and son. Through the accurate identification and identification of the above-mentioned characteristic sequences, the genetic material of the offspring organism can be effectively extracted without a reference genome and a high-quality public database. Therefore, in this step, based on the sequences shared by these parents or unique to the father or the mother, the parent-common characteristic fragment set is defined as a set consisting of sequence fragments shared by the father and mother of the organism, the father-specific characteristic fragment set is defined as a set consisting of sequence fragments unique to the father of the organism, and the mother-specific characteristic fragment set is defined as a set consisting of sequence fragments unique to the mother of the organism.

[0064] In some embodiments, a set operation is performed based on the paternal Kmer sequence set and the maternal Kmer sequence set to obtain a characteristic fragment set including a parental common characteristic fragment set, a paternal unique characteristic fragment set, and a maternal unique characteristic fragment set.

[0065] Kmer (K-mer) sequence refers to a sequence of length k intercepted from any nucleic acid sequence (read) in the sequencing data. For a nucleic acid sequence (read) with a read length of n, starting from the first base and sliding with a step length of one base, (n-k+1) Kmer sequences of length k can be obtained.

[0066] Wherein, set refers to the whole of several elements that are determined and can be distinguished. In the present application, element can be the read length obtained by the above-mentioned sequencing or a partial sequence in the read length, such as a Kmer sequence. Therefore, the paternal Kmer sequence set refers to all Kmer sequences obtained by all nucleic acid sequences measured in the sequencing data of the paternal parent, and the maternal Kmer sequence set refers to all Kmer sequences obtained by all nucleic acid sequences measured in the sequencing data of the maternal parent.

[0067] Among them, set operation refers to forming a new set through operation of two or more sets, and the specific operation methods include at least one of the union, intersection, difference, complement, symmetric difference, etc. In this scheme, the set operation method is mainly used to obtain the common feature fragment set of parents, the unique feature fragment set of the father, and the unique feature fragment set of the mother, so its set operation method includes at least one of the intersection and difference. Taking set A and set B as an example, the intersection A∩B is a new set composed of the common elements of set A and set B; the difference set AB is a new set composed of elements belonging to A but not to B; the difference set BA is a new set composed of elements belonging to B but not to A. For example, based on the paternal Kmer sequence set and the maternal Kmer sequence set, the common feature fragment set of parents is obtained by the set operation of intersection, and the paternal unique feature fragment set and the maternal unique feature fragment set are obtained by the set operation of difference.

[0068] In some embodiments, the length of the Kmer sequence used to construct the characteristic fragment set is 10 to 100 bases, for example, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48 , 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100 bases. In some embodiments, the length of the Kmer sequence used to construct the characteristic fragment set is 10-90, 10-80, 10-70, 10-60, 10-50, 10-40, 10-30 bases.

[0069] In some embodiments, the paternal Kmer sequence set and the maternal Kmer sequence set are the paternal Kmer sequence set and the maternal Kmer sequence set after deleting the sequencing error sequence and the repetitive sequence. Due to the contingency and randomness of sequencing errors, the number of occurrences of the Kmer sequence containing a specific error base is much lower than the Kmer sequence obtained by normal sequencing; and the repetitive sequence is the opposite, the number of occurrences of the corresponding Kmer sequence is much higher than the Kmer sequence obtained by normal sequencing, so the Kmer sequences in the paternal Kmer sequence set and the maternal Kmer sequence set are respectively subjected to frequency statistics, and the sequencing error peaks with a frequency lower than the ultra-low frequency threshold and the repetitive sequence peaks with a frequency higher than the ultra-high frequency threshold are deleted. In some embodiments, the frequency thresholds of ultra-low frequency and ultra-high frequency are automatically obtained by the Kmer frequency distribution diagram. For example, according to the Kmer frequency distribution diagram, the ultra-low frequency domain value is defined as the frequency corresponding to the first trough (the dividing line between the sequencing error peak and the main peak); the ultra-high frequency domain value is defined as the frequency of multiples (such as 5 times) of the frequency value corresponding to the highest peak (sequencing main peak).

[0070] In some embodiments, considering the sequencing cost, the paternal sequencing data and the maternal sequencing data use relatively inexpensive second-generation whole genome sequencing data (NGS) instead of relatively expensive single molecule sequencing data.

[0071] In actual production, due to the use of different long-read sequencing technology platforms, the difficulty of obtaining characteristic sequence information is significantly different. The main reason is that the difference in single-base error rate and read length distribution of different sequencing technologies brings great obstacles and challenges to the accurate acquisition of characteristic sequence information. Therefore, according to the difference in sequencing error rate and length, almost all single-molecule long-read sequencing technology platforms on the market can be classified: that is, one is a long-read sequencing technology with a shorter length but a lower error rate, such as the high-fidelity long-read sequencing data (PacBio HiFi) of the Pacific Biosciences (Pacific Biosciences) technology platform; one is a long-read sequencing technology with a longer length but a higher error rate, such as the Continuous Long Reads (PacBio CLR) data of the Pacific Biosciences (Pacific Biosciences) technology platform, the ONT sequencing data of Oxford Nanopore Technologies, and the long-read sequencing data of the Cyclone WT platform of the domestic sequencing platform of Shenzhen BGI Life Sciences Research Institute. For this reason, two different methods are proposed in the embodiment of the present application to obtain the genetic characteristic sequence information of an organism: one is based on the Kmer fragment unit, and the other is based on the Strobemer fragment unit.

[0072] In some embodiments, for the sequencing data of an organism obtained by a long read sequencing technology with a low sequencing error rate, such as PacBio HiFi sequencing data, Kmer sequence units are used to efficiently and accurately identify genomic features (such as genomic repeats and heterozygous sequences), and the information of family genetic feature sequences is quickly obtained. Therefore, the feature fragment is a Kmer sequence, and the feature fragment set is a Kmer sequence set. Set operations are directly performed with the paternal Kmer sequence set and the maternal Kmer sequence set to obtain the parental Kmer sequence set, the paternal Kmer sequence set, and the maternal Kmer sequence set as the parental Kmer sequence set, respectively. The parental Kmer sequence set is recorded as set P, and the maternal Kmer sequence set is recorded as set M. Set P is used to take a difference set to set M, that is, the set of Kmer sequences contained in set P but not in set M, and finally the paternal Kmer sequence set POK is obtained. Use set M to take the difference set of set P, that is, the set of Kmer sequences contained in set M but not in set P, and finally get the mother-specific Kmer sequence set MOK. The intersection of set P and set M is the parent-shared Kmer sequence set SK.

[0073] Although some long-read sequencing technologies can obtain ultra-long fragment sequencing results, they cannot always maintain a high level of sequencing accuracy, resulting in a single-base error rate that is much higher than that of traditional second-generation sequencing technologies. When using Kmer sequence units to directly analyze the long reads of these offspring, more low-frequency Kmer sequence units are generated. The presence of these Kmer sequences will cause distortion of genomic characteristic peaks, and effective Kmer sequences carrying paternal and maternal specific heterozygous site information will also be mistakenly classified as low-frequency areas due to the presence of erroneous bases, and will be mistakenly identified as erroneous peaks and removed. Therefore, accurate Kmer matching technology is obviously unable to analyze such long-read sequences with high sequencing error rates, which brings difficulties to analyzing family genome characteristics and identifying offspring long reads. To this end, in some embodiments, for the sequencing data of an organism obtained by a long read sequencing technology with a longer length but a higher error rate, such as PacBio CLR sequencing data, ONT sequencing data, and Cyclone WT sequencing data, the set operation also includes using the Strobemer method to convert the paternal Kmer sequence set and the maternal Kmer sequence set into the paternal Strobemer sequence set and the maternal Strobemer sequence set, and innovatively applying the Strobemer sequence unit with fault tolerance to extract the family genetic characteristic sequence information from the long read sequence of the offspring organism containing more sequencing errors. Therefore, in these embodiments, the characteristic fragment is a Strobemer sequence, and the characteristic fragment set is a Strobemer sequence set. The parental common characteristic fragment set generated by the set operation is the parental common Strobemer sequence set, the paternal unique characteristic fragment set is the paternal unique Strobemer sequence set, and the maternal unique characteristic fragment set is the maternal unique Strobemer sequence set. Similar to the Kmer sequence method, the paternal Strobemer sequence set is recorded as set P, and the maternal Strobemer sequence set is recorded as set M. Use set P to take the difference set of set M to obtain the father's unique Strobemer sequence set POK. Use set M to take the difference set of set P to obtain the mother's unique Strobemer sequence set MOK. The intersection of set P and set M is the parental common Strobemer sequence set SK.

[0074] The Strobemer method is a method based on the combination of Kmer error-tolerant compression tag sequence information. The general method of Strobemer extraction is: for any sequence of sufficient length (sufficient to extract at least one Strobemer), based on variable size (w min ,w max ) sliding window, extracting an l-mer from each window, connecting n consecutive l-mers head to tail, and the obtained n*l-mer is used as the representative strobemer of the sequence covered by this continuous sliding window.

[0075] In some embodiments, l is 10-100, for example, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100; n is 2-10, for example, 2, 3, 4, 5, 6, 7, 8, 9, 10; w min is 5 to 20, for example, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20; w max is 10 to 50, for example, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 25, 30, 35, 40, 45, 50; and w min <w max .

[0076] In some embodiments, the Strobemer method includes at least one of minstrobes, hybridstrobes, randstrobes, mixedstrobes, altstrobes, multistrobes, etc. The specific operating methods and principles of these different Strobemer methods can at least refer to Sahlin, K. (2021). Effective sequence similarity detection with strobemers. Genome Research, gr-275648.; and Maier, BD, & Sahlin, K. (2023). Entropy predicts sensitivity of pseudo-random seeds. Genome Research, gr-277645.

[0077] Among them, clustering refers to an unsupervised learning method that divides multiple objects into different groups (or clusters) through a static classification method so that member objects in the same group have similar attributes. Clustering includes clustering methods based on partitioning (such as k-means, k-means++, bi-kmeans, FCM, etc.), hierarchy (such as single link, full link, even link, etc.), density (DBSCAN, OPTICS), grid (such as STING, CLIQUE), model (CLASSI, Cheeseman, AutoClass), graph, etc. In some implementation methods, clustering uses a Bayesian-Gaussian mixture model.

[0078] The general process of clustering is divided into data preparation, feature selection, feature extraction, clustering, and result evaluation. In the clustering method of the Bayesian-Gaussian mixture model, the clustering process includes dimensionality reduction. For example, through linear dimensionality reduction principal component analysis, the original feature variables are transformed to form independent principal components, objectively and truly capturing the elements with higher contribution; according to the data density distribution characteristics after dimensionality reduction, the Gaussian probability density function (normal distribution curve) is used to quantify the principal components, and the Bayesian-Gaussian mixture model is applied for clustering. The algorithm uses the initial unlabeled data to solve the prior probability of the parameters of each Gaussian model; then the maximum likelihood function is used to iteratively solve, and the maximum a posteriori probability class is found through the Bayesian method. The algorithm can effectively reduce the interference of abnormal data noise points caused by high error rates, reduce dependence on initial parameters, and benefit the development of unknown species projects.

[0079] Among them, species characteristics refer to sequence-related parameters with species specificity. In the sequencing process, in order to obtain effective genetic information, the read length with a depth coverage of more than ten or dozens of times the species genome is usually measured. Since the obtained offspring read length sequence generally has a large amount of data, and the read length sequence obtained by the long read length sequencing technology is not fixed in length, and is accompanied by problems such as high error rate, it is more difficult to obtain the offspring sequence characteristics. The combined species difference characteristics can effectively improve the anti-interference ability of family model pollution filtering, avoid information misleading caused by high error rate, and improve the efficiency of offspring identification that only relies on family characteristic sequences. In order to make full use of the unique GC non-preference and long read length advantages of the third-generation data, in some embodiments, at least one of the GC content and oligonucleotide markers is used as a non-family species characteristic.

[0080] Among them, GC content refers to the percentage of the sum of the number of guanine and cytosine in a nucleotide sequence to the total number of bases in the sequence. GC content is one of the better genetic reference standards for species classification and identification. For example, the base composition of microbial DNA is species-specific, that is, the (G+C) value of microbial DNA is extremely stable in cells and is not affected by factors other than bacterial age and mutation factors. The GC content of DNA molecules of different species is relatively stable. Therefore, in some embodiments, GC content can be used as a species feature on which the residual sequence clustering depends, recorded as a statistic of feature x4. It is understandable that during the operation, for global statistical information such as GC content, the read length features are combed, the functional relationship between species discrimination and read length is calculated, and the minimum length threshold of tolerable erroneous bases is determined. It should be noted that, affected by the GC preference of different sequencing technologies, GC content may introduce the characteristics of sequencing technology to a certain extent. Therefore, in some further embodiments, species characteristics need to be comprehensively corrected by introducing other information. Therefore, species characteristics also include oligonucleotide markers.

[0081] The oligonucleotide identifier of the read length refers to a Kmer sequence of a specific order. In some embodiments, the oligonucleotide identifier of the read length is represented by the frequency of occurrence of the Kmer sequence of a specific order. For a Kmer sequence of length K, there are 4 K (K>=2), for example, when K is 2, the types of oligonucleotide identifiers are 4 to the power of 2, and 16 types of fragment units are obtained; when K is 3, the types of oligonucleotide identifiers are 4 to the power of 3, and 64 types of fragment units are obtained. Therefore, for a read length sequence of length L, the frequency of occurrence of 16 oligonucleotide identifiers or 64 oligonucleotide identifiers in its L-K+1 Kmer sequence can be obtained. In some embodiments, the length of the oligonucleotide identifier is 2 to 10, for example, 2nt, 3nt, 4nt, 5nt, 6nt, 7nt, 8nt, 9nt, 10nt. It is understandable that in practice, if the K value is too small, the differences between species cannot be effectively enumerated; if the K value is too large, it will lead to a large consumption of computing resources, resulting in low efficiency in species difference identification. Therefore, based on the particularity of the genetic codon and the rationality of computing resource allocation, the fragment unit length K is generally used as 3. In addition, due to the double helix structure of DNA and the characteristics of base pairing, in actual sequencing, it is impossible to predict whether the sequenced strand of the two complementary strands is the template strand. In actual statistics, the two complementary fragment units are combined for statistical calculation. In some embodiments, the frequency of all fragment units is counted for each read length and recorded as the features x5, x6, ..., x(4 K / 2+4). For local base information such as oligonucleotide markers, the species differentiation efficiency and accuracy of multiple short oligonucleotide lengths are tested according to the species genome size to reduce the sensitivity to incorrect bases. Therefore, in some other embodiments, for oligonucleotide markers, the frequency of fragment units of some Kmer sequences can be selected for statistics or the fragment units of Kmer sequences of different lengths can be comprehensively counted or a combination of the two.

[0082] By statistically analyzing the distribution of each differential feature in the host and contaminant species, the incompleteness and complementarity of its contribution to contamination filtering have been found: the family model is affected by the read error rate and length distribution, and it is impossible to complete the labeling of all host data; the nuclear extraction and sterilization features are limited by technology, and only part of the prokaryotic contamination is labeled; GC and oligonucleotide markers are also susceptible to the characteristics of the third generation of data. Therefore, it is necessary to construct a high-dimensional matrix of multiple differential features to utilize the complementary relationship between features and jointly participate in the distinction of read lengths in the clustering hyperspace. Therefore, in some of the embodiments, for the residual sequences, clustering is further performed through the characteristic GC content and the characteristic joint system of oligonucleotide markers to further improve the sensitivity of host sequence recognition. In other embodiments, for the residual sequences, the features selected for clustering can include the content of the characteristic fragment set in addition to GC and oligonucleotide markers, such as at least one of the content of the parental common characteristic fragment set, the content of the paternal-specific characteristic fragment set, and the content of the maternal-specific characteristic fragment set.

[0083] In view of the different distribution characteristics of the three generations of read lengths, all features were normalized. At the same time, the residual Kmer marker features after centrifugal cell nucleus extraction and sterilization were used to further integrate the application of bioinformatics and experimental technology features in the pollution filtration system.

[0084] S120: Marking a target read group having the high-confidence read sequence from a plurality of read groups, wherein the high-confidence read sequence is a sequence whose content of the characteristic fragment set is above a threshold.

[0085] Among them, the content refers to the ratio of the characteristic fragments contained in the read sequence of the sequencing data to the total number of characteristic fragments in the characteristic fragment set. For example, the content of the paternal-specific characteristic fragment set of a certain read sequence in the sequencing data refers to the ratio of the paternal-specific characteristic fragments contained in the read sequence to the total number of characteristic fragments in the paternal-specific characteristic fragment set, recorded as feature x1; the content of the maternal-specific characteristic fragment set refers to the ratio of the maternal-specific characteristic fragments contained in a certain read sequence of the sequencing data to the total number of characteristic fragments in the maternal-specific characteristic fragment set, recorded as feature x2; the content of the parental common characteristic fragment set refers to the ratio of the parental common characteristic fragments contained in a certain read sequence of the sequencing data to the total number of characteristic fragments in the parental common characteristic fragment set, recorded as feature x3.

[0086] In some embodiments, the content of the characteristic fragment set is above the threshold value, which means that the content of the common characteristic fragment set of the parents, the content of the paternal-specific characteristic fragment set and the content of the maternal-specific characteristic fragment set are above the threshold value. The sequence in which the content of the characteristic fragment set is above the threshold value refers to the read sequence in which the content of the paternal-specific characteristic fragment set is above the first threshold value, the content of the maternal-specific characteristic fragment set is above the second threshold value, and the content of the common characteristic fragment set of the parents is above the third threshold value. The first threshold value, the second threshold value, and the third threshold value may be the same or different. In some embodiments, it may be required that the sum of the content of the paternal-specific characteristic fragment set and the content of the maternal-specific characteristic fragment set is above the fourth threshold value as a further limitation, or as a substitute for the content of the paternal-specific characteristic fragment set being above the second threshold value and the content of the maternal-specific characteristic fragment set being above the third threshold value. In some embodiments, it may be required that the sum of the content of the common characteristic fragment set of the parents, the content of the paternal-specific characteristic fragment set and the content of the maternal-specific characteristic fragment set be above the fifth threshold value. In other embodiments, the content of the characteristic fragment set is above the threshold value may also mean that at least one of the content of the parental common characteristic fragment set, the content of the paternal unique characteristic fragment set and the content of the maternal unique characteristic fragment set is above the threshold value.

[0087] In some embodiments, the sequence in which the content of the parental common characteristic fragment set, the content of the paternal specific characteristic fragment set and the content of the maternal specific characteristic fragment set are above the threshold refers to the read sequence in which the content of the paternal specific characteristic fragment set is greater than the first threshold, the content of the maternal specific characteristic fragment set is greater than the second threshold, and the content of the parental common characteristic fragment set is greater than the third threshold. In some embodiments, it can also be required that the sum of the content of the paternal specific characteristic fragment set and the content of the maternal specific characteristic fragment set is greater than the fourth threshold as a further limitation, or as a substitute for the content of the paternal specific characteristic fragment set being greater than the second threshold and the content of the maternal specific characteristic fragment set being greater than the third threshold. In some embodiments, it can also be required that the sum of the content of the parental common characteristic fragment set, the content of the paternal specific characteristic fragment set and the content of the maternal specific characteristic fragment set is greater than the fifth threshold.

[0088] In some embodiments, the first threshold, the second threshold, the third threshold, the fourth threshold, and the fifth threshold are independently 0.01-1%, for example, 0.01%, 0.02%, 0.03%, 0.04%, 0.05%, 0.06%, 0.07%, 0.08%, 0.09%, 0.1%, 0.2%, 0.3%, 0.4%, 0.5%, 0.6%, 0.7%, 0.8%, 0.9%, 1%. In some embodiments, the first threshold is independently 0.05-0.5%, the second threshold is independently 0.05-0.5%, and the third threshold is independently 0.2-0.8%. In some embodiments, the first threshold is independently 0.05-0.2%, the second threshold is independently 0.05-0.2%, and the third threshold is independently 0.4-0.6%.

[0089] In some embodiments, it also includes determining whether the read sequence in the sequencing data is a paternal-specific read sequence or a maternal-specific read sequence according to the content of the paternal-specific characteristic fragment set and the content of the maternal-specific characteristic fragment set. For example, if the content of the paternal-specific characteristic fragment set is above or greater than the sixth threshold, the read sequence is determined to be a paternal-specific read sequence; if the content of the maternal-specific characteristic fragment set is above or greater than the seventh threshold, the read sequence is determined to be a maternal-specific read sequence. In this way, the read sequence of the offspring is classified as paternal or maternal-specific through features x1 and x2, which can lay the foundation for the next step of fully typing the haploid genome assembly of the host species.

[0090] Among them, the paternal-specific characteristic fragment set POK, the maternal-specific characteristic fragment set MOK, and the parent-shared characteristic fragment set SK extracted using the characteristic fragment unit information (Kmers or Strobemer) are used as family marker sequence information, and the family information density in the read sequence of the offspring is calibrated using the marker sequences of the three sets to obtain the three-dimensional statistical data of each offspring read sequence, which are recorded as the statistics of features x1, x2, and x3 respectively. In some embodiments of the present application, the provided sequencing data filtering method is to filter based on the family model, extracting the offspring that contains a considerable number of paternal or maternal-specific and parent-shared fragment units (Kmer sequence units or Strobemer sequence units), that is, the features x1, x2, and x3 must simultaneously exceed their respective set thresholds, and can be marked as high-confidence read sequences belonging to the organism itself in the organism sequencing data.

[0091] After clustering is completed to obtain multiple groups, the target group is identified based on the high-confidence read sequence. Since the core of the algorithm adopts an unsupervised clustering algorithm that does not introduce any prior knowledge, how to host-label each cluster group of clustering results is one of the main issues. In the model, features x1 to x3 should only be enriched in the read group belonging to the organism itself. Based on this, this is applied in the present invention to further identify the sequence cluster groups of organisms for all unsupervised clustering results.

[0092] In some embodiments of the present application, the target read group is a read group whose high-confidence read sequence is greater than or equal to a critical value. For example, the proportion of high-confidence read sequences in each read group or the content density distribution of the characteristic fragment set of each read sequence is counted, and the read groups greater than or equal to the proportion threshold or density distribution threshold are screened out and marked as the target read group of the organism to be selected. In addition, the proportion of high-confidence read sequences can be sorted from large to small, and the first 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 sequence groups are selected as the target read groups marked this time.

[0093] S130: Repeat the cycle of S110 and S120 to screen out read length sequences whose occurrence frequency in all target read length groups is higher than a set value.

[0094] Due to the randomness of the unsupervised clustering algorithm and the complexity of the impact of sequence characteristics on sequencing errors, it is not possible to ensure that some host long reads are enriched in the same result group each time, affecting the integrity and purity of the host species cluster group, and thus affecting the effect and efficiency of the final sequence decontamination. Therefore, multiple cycles of unsupervised clustering can be used, and then the target read length group marked by each result can be integrated.

[0095] In some embodiments, the number of cycles is 5 to 20, for example, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20. In some embodiments, the set value is 60 to 99%, for example, 60%, 61%, 62%, 63%, 64%, 65%, 66%, 67%, 68%, 69%, 70%, 71%, 72%, 73%, 74%, 75%, 76%, 77%, 78%, 79%, 80%, 81%, 82%, 83%, 84%, 85%, 86%, 87%, 88%, 89%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99%.

[0096] In some specific embodiments of the present invention, all read sequences are clustered using a Bayesian-Gaussian mixture model based on the features of features x1 to x36, the sequence groups obtained by clustering are sorted according to the content of the feature fragment set containing features x1 to x3, and the first 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 sequence groups are selected as the target read groups for this time; the cycle is repeated 5 to 20 times, and finally a set of a specified number of target read groups can be obtained; the frequency of each read sequence appearing in the set of all target read groups is counted, and the read sequences above the set value are integrated into the final total reads of the organism.

[0097] In the implementation of the present invention, the reads can be classified according to the labels of each data set through the aforementioned steps, and the data collection of the biological host and the symbiotic parasitic microorganisms can be realized by identifying the number of each read group and the Kmer peak of each read group, and then targeted single species or multi-species mixed assembly software can be selected for assembly respectively.

[0098] Therefore, a second aspect of the present invention provides a method for assembling sequencing data, which comprises filtering the sequencing data according to the aforementioned method for filtering sequencing data of an organism before assembly.

[0099] In some embodiments, the enrichment and assembly of the read sequences of the organism includes assembling all the read sequences of the final organism obtained by filtration through a genome assembly tool. In some embodiments, the assembly tool is selected from single-molecule long-read benchmark genome assembly tools such as Canu and Flye. Specifically, all the read sequences of the final organism are extracted and stripped from the original sequencing data, and then mixed into an input file. The input file is further assembled using single-molecule long-read benchmark genome assembly tools such as Canu and Flye.

[0100] In some embodiments, if the sequencing depth, read length, host species heterozygosity, etc. of the read sequence meet the requirements, the extracted read sequence of the organism can be further decomposed into two haplotype read sequence groups, namely, the paternal-specific haplotype read sequence group and the maternal-specific haplotype read sequence group, through the paternal-specific characteristic fragment set and the maternal-specific characteristic fragment set. Then, two fully typed haplotype genomes are produced using genome assembly tools.

[0101] In some embodiments, after extracting all read sequences of an organism from its sequencing data, the remaining sequences, such as the long reads of contaminants such as symbionts, still have research value for analyzing the host's living environment, symbiotic and parasitic relationships, etc. Therefore, the separated remaining sequences can be further explored for their biological components. The remaining sequences are extracted from the original data, and then analyses such as single-molecule long read assembly and variation can be performed.

[0102] According to a third aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores computer-executable instructions for causing a computer to execute the aforementioned filtering method or assembling method.

[0103] A fourth aspect of the present invention provides an electronic device, comprising a processor and a memory, wherein the memory stores a computer program executable on the processor, and the processor implements the aforementioned filtering method or assembling method when running the computer program.

[0104] The memory, as a non-transient computer-readable storage medium, can be used to store non-transient software programs and non-transient computer executable programs, such as the aforementioned filtering method or assembly method described in the embodiments of the present application. The processor implements filtering or assembly of sequencing data by running the non-transient software programs and instructions stored in the memory.

[0105] The memory may include a program storage area and a data storage area, wherein the program storage area may store an operating system and an application required for at least one function; and the data storage area may store and execute the above-mentioned programs. In addition, the memory may include a high-speed random access memory and may also include a non-volatile memory, such as at least one disk storage device, a flash memory device, or other non-volatile solid-state storage devices.

[0106] In some embodiments, the memory may include a memory remotely located relative to the processor, and the remote memory may be connected to the processor via a network. Examples of the above network include, but are not limited to, the Internet, an intranet, a local area network, a mobile communication network, and combinations thereof.

[0107] The non-transitory software program and instructions required to implement the above method are stored in the memory, and when executed by one or more processors, the above filtering method or assembling method is performed.

[0108] A fifth aspect of the present invention provides a sequencing data filtering system, comprising:

[0109] A clustering module, wherein the clustering module is used to cluster the sequencing data of the organism according to the content and species characteristics of the characteristic fragment set to obtain a plurality of read length groups, wherein the characteristic fragment set includes at least one of a parental common characteristic fragment set, a paternal unique characteristic fragment set, and a maternal unique characteristic fragment set;

[0110] A marking module, the marking module is used to mark a target read length group having the high-confidence read length sequence from the plurality of read length groups, the high-confidence read length sequence being a sequence whose content of the characteristic fragment set is above a threshold value;

[0111] A screening module is used to screen out read sequences whose occurrence frequency in all target read length groups is higher than a set value after the clustering module and the marking module are repeatedly cycled.

[0112] In some embodiments, the characteristic fragment set is obtained based on the set operation of the father Kmer sequence set and the mother Kmer sequence set.

[0113] In some embodiments, the set operation includes taking the intersection and difference of two sets.

[0114] In some embodiments, the characteristic fragment set obtained based on the set operation of the paternal Kmer sequence set and the maternal Kmer sequence set includes:

[0115] Based on the sequencing data of the paternal and maternal copies of the organism, a paternal Kmer sequence set and a maternal Kmer sequence set are obtained;

[0116] A set operation is performed based on the paternal Kmer sequence set and the maternal Kmer sequence set to obtain a characteristic fragment set.

[0117] In some of the embodiments, the intersection of the paternal Kmer sequence set and the maternal Kmer sequence set is taken to obtain the parental common characteristic fragment set; the difference between the paternal Kmer sequence set and the maternal Kmer sequence set is taken to obtain the paternal unique characteristic fragment set and the maternal unique characteristic fragment set.

[0118] In some embodiments, the characteristic fragment set is a Kmer sequence set.

[0119] In some embodiments, the length of the Kmer sequence in the Kmer sequence collection is 10 to 100 bases.

[0120] In some embodiments, the set operation also includes using the Strobemer method to convert the paternal Kmer sequence set and the maternal Kmer sequence set into the paternal Strobemer sequence set and the maternal Strobemer sequence set; the characteristic fragment set is the Strobemer sequence set.

[0121] In some embodiments, the set operation includes taking the intersection and difference of the father Strobemer sequence set and the mother Strobemer sequence set.

[0122] In some embodiments, the intersection of the paternal Strobemer sequence set and the maternal Strobemer sequence set is taken to obtain the parental common characteristic fragment set; the difference between the paternal Strobemer sequence set and the maternal Strobemer sequence set is taken to obtain the paternal unique characteristic fragment set and the maternal unique characteristic fragment set.

[0123] In some embodiments, the strobemer sequence in the strobemer sequence set is based on a sliding window of variable size, each window extracts an l-mer, and connects n consecutive l-mers head to tail to obtain n*l-mer; wherein l is 10 to 100, n is 2 to 10, and the variable size has a minimum value w min and the maximum value w max , w min 5~20,w max is 10 to 50, and w min <w max .

[0124] In some embodiments, the strobemer method includes at least one of minstrobes, hybridstrobes, randstrobes, mixedstrobes, altstrobes, and multistrobes.

[0125] In some embodiments, the sequencing data of the paternal and maternal copies of the organism are second-generation sequencing data.

[0126] In some of the embodiments, the sequencing data of the organism is third-generation sequencing data.

[0127] In some embodiments, the third-generation sequencing data is any one of PacBio HiFi, PacBio CLR, ONT, and CycloneWT sequencing data.

[0128] In some embodiments, the sequence whose characteristic fragment set content is above the threshold is a sequence whose father-specific characteristic fragment set content is above the first threshold, whose mother-specific characteristic fragment set content is above the second threshold, and whose parent-common characteristic fragment set content is above the third threshold.

[0129] In some embodiments, the first threshold, the second threshold and the third threshold are independently 0.01-1%.

[0130] In some embodiments, the first threshold is independently 0.05-0.5%, the second threshold is independently 0.05-0.5%, and the third threshold is independently 0.2-0.8%.

[0131] In some embodiments, the first threshold is independently 0.05-0.2%, the second threshold is independently 0.05-0.2%, and the third threshold is independently 0.4-0.6%.

[0132] In some embodiments, the species characteristic comprises at least one of GC content and frequency of oligonucleotide markers.

[0133] In some embodiments, the length of the oligonucleotide tag is 2 to 10.

[0134] In some embodiments, the length of the oligonucleotide tag is 2 to 4.

[0135] In some embodiments, the length of the oligonucleotide tag is 3.

[0136] In some embodiments, the frequencies of oligonucleotide markers are combined with the frequencies of complementary oligonucleotides.

[0137] In some embodiments, the number of cycles is 5 to 20 times.

[0138] In some embodiments, the set value is 60-99%.

[0139] In some embodiments, the assembly tool is selected from Canu, Flye.

[0140] In some of the embodiments, the method further includes classifying the read sequences into paternal-specific read sequences and maternal-specific read sequences according to the paternal-specific characteristic fragment set and the maternal-specific characteristic fragment set, and assembling them separately.

[0141] The device or system implementation described above is only illustrative, and the modules described as separate components may or may not be physically separated, that is, they may be located in one place or distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0142] It is understood that all or some of the steps disclosed above, etc. can be implemented as software, firmware, hardware and appropriate combinations thereof. Some physical components or all physical components can be implemented as software executed by a processor, such as a central processing unit, a digital signal processor or a microprocessor, or implemented as hardware, or implemented as an integrated circuit, such as an application-specific integrated circuit. Such software can be distributed on a computer-readable medium, and the computer-readable medium can include a computer storage medium (or a non-transitory medium) and a communication medium (or a temporary medium). It is understood that computer storage media include volatile and non-volatile, removable and non-removable media implemented in any method or technology for storing information (such as computer-readable instructions, data structures, program modules or other data). Computer storage media include, but are not limited to, RAM, ROM, EEPROM, flash memory or other memory technology, CD-ROM, digital versatile disk (DVD) or other optical disk storage, magnetic cassettes, magnetic tapes, disk storage or other magnetic storage devices, or any other medium that can be used to store desired information and can be accessed by a computer.

[0143] Additionally, it will be appreciated that communication media typically embodies computer readable instructions, data structures, program modules, or other data in a modulated data signal such as a carrier wave or other transport mechanism, and may include any information delivery media.

[0144] The present invention is further described in detail below through specific examples.

[0145] It should be understood that these embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention.

[0146] The experimental methods in the following examples without specifying specific conditions are usually carried out under conventional conditions or under conditions recommended by the manufacturers. The materials and reagents used in the examples are commercially available unless otherwise specified.

[0147] Example 1: Long-read sequencing simulation data of artificial host-commensal mixing of three lineages

[0148] This embodiment uses samples with family composition in the 1,000 Genomes to construct a human-intestinal microbial symbiont model. The family information is selected as the paternal number HG003 and the maternal number HG004 as the parental source of the family sample; the offspring number is HG002 as the host organism to construct artificially simulated symbiont data. The NGS data volume of the father is 75,327,000 paired-end (PE) sequences, and the NGS data volume of the mother is 75,396,412 paired-end sequences as the parent data set. The reference genome information of the offspring sample HG002, human intestinal microorganisms, random contamination sequences, and the reference genome of the mouse of a closely related species are used as the mixed sequences of the offspring DNA sequencing as the simulation objects, and the long read data simulation tool is used to simulate the three most mainstream long read sequencing technologies on the market (PacBio HiFi, PacBio CLR, ONT), and finally different types of data sets are generated. These data features include sequencing platform, average read length, average error rate, etc. The specific data simulation features are as follows:

[0149] PacBio HiFi simulation data set: human chromosome 19 in HG002 is used as the host DNA source, two bacteria from the unified human gastrointestinal genome (UHGG) collection are used as interphylum symbionts, and another 8 UHGG bacteria are used as random contaminants. Mouse chromosome 7 is used as a symbiont that shares more similar sequences with the host genome to test the sensitivity to the similarity of host and symbiotic genome sequences. PBSIM2 simulates PacBio HiFi long reads based on the reference genome, using the pre-trained read length and single-base error rate distribution for the PacBio HiFi sequencing platform provided by PBSIM2, with an average read length of 10 kilobases (kbp) and an average error rate of 1%. A total of 289,446 entries.

[0150] PacBio CLR simulation data: The components are the same as those of PacBio HiFi simulation data. The pre-trained read length and single-base error rate distribution for the PacBio CLR sequencing platform provided by PBSIM2 are used, with an average read length of 10kbp and an average error rate of 5%. A total of 285,181 records.

[0151] ONT simulation data: The composition is the same as that of PacBio HiFi simulation data. The pre-trained read length and single-base error rate distribution for the ONT sequencing platform provided by PBSIM2 are used, with an average read length of 30kbp and an average error rate of 15%. A total of 96,575 records.

[0152] For all simulations, the default error profile provided by PBSIM2 was used (PacBio substitution: insertion: deletion = 6:50:54, ONT 23:31:46). 50× host data were simulated using chromosome 19 of the HG002 reference genome (GCA_011064465.1), and the quality of the simulated data is shown in Table 1. 50× mouse chromosome 7 data were simulated based on the GRCm39 reference genome (GCA_000001635.9). Other bacteria were simulated based on the UHGG reference genome with various coverages. See Table 1 for details.

[0153] Table 1. Simulated data component information table of PacBio HiFi, PacBio CLR and ONT sequencing modes

[0154]

[0155] 1. Extraction of family genetic characteristic sequence information (characteristic fragment set)

[0156] 1. Extraction of family genetic feature sequence information based on K-mer units

[0157] When constructing the genetic characteristics information of the family, the K-mer generation tool (meryl tool v1.3) was first used to K-merize the NGS data of the paternal HG003 and maternal HG0004. Considering the size of the human genome, the length and error rate distribution of long reads, and the heterozygosity of the human genome, K=21-mer was adopted. The meryl tool was used to perform set operations on the rapidly generated 21-mers to generate the paternal-specific set POK, the maternal-specific set MOK, and the 21-mer set SK shared by the two. The specific process is as follows: the meryl tool is used to extract the 21-mer fragment unit of each read length in the paternal and maternal sequencing data, that is, firstly, based on the sequence fragment unit with a fixed length of 21 bases, each read length of the paternal and maternal is started from the first base, and the sliding window form (step size is 1) is used to extract the fragment units of length 21 in sequence, and the 21-mer fragment units of all read lengths are obtained in sequence. Secondly, the 21-mer sequence frequencies of the father and mother are counted respectively, and the sequencing error peaks and repeated sequence peaks are deleted using the frequency statistical distribution, that is, the ultra-low frequency and ultra-high frequency 21-mer sequences are deleted. At this time, the father's 21-mer set is recorded as set P, and the maternal 21-mer set is recorded as set M. Then, in order to obtain the father's unique 21-mer set, set P is used to take the difference set of set M, and finally the father's unique set POK is obtained. In order to obtain the mother's unique 21-mer set, set M is used to take the difference set of set P, and finally the mother's unique set MOK is obtained. Finally, the intersection of set P and set M is taken, which is the 21-mer set shared by the father and mother, recorded as SK. Finally, the father's and mother's unique 21-mer sets POK, MOK and the total 21-mer set SK are obtained.

[0158] Finally, the characteristic sequence information was obtained to establish a family characteristic set with a K-mer fragment unit of 21bp in length, with a total of 8426128 21-mers in the paternal-specific set POK, 8695115 21-mers in the maternal-specific set MOK, and 163185258 21-mers in the parent-common set SK. The total number of 21-mers marked in the family was 180306501.

[0159] 2. Extraction of characteristic information of pedigree genetic sequences based on strobemer unit fragments

[0160] Use Strobemer fragment units to extract family characteristic sequence information. After a certain amount of data pre-testing, this scheme chooses to use random Strobes with n=2, w_min=15, and w_max=30 as a practical technology. First, obtain a 40-mer K-mer feature database. Secondly, convert the 40-mers one by one into corresponding Strobemer representative sequences. At this time, the paternal Strobemer set is recorded as set P, and the maternal Strobemer set is recorded as set M. Then, in order to obtain the paternal-specific Strobemer set, use set P to take the difference set of set M, and finally obtain the paternal-specific set POK. In order to obtain the maternal-specific Strobemer set, use set M to take the difference set of set P, and finally obtain the maternal-specific set MOK. Finally, the intersection of set P and set M is the paternal and maternal common Strobemer set, recorded as SK. Finally, the parent-specific Strobemer sets POK, MOK and the common Strobemer set SK are obtained.

[0161] Finally, a family characteristic set of Strobemer fragment units carrying characteristic sequence information was obtained, with a total of 27,794,812 Strobemers in the paternal-specific set POK, 26,295,868 Strobemers in the maternal-specific set MOK, and 325,180,216 Strobemers in the parent-common set SK. The total number of Strobemers marked in the family was 379,270,896.

[0162] In the extraction of family characteristics, the paternal-specific set POK, the maternal-specific set MOK, and the parent-common set SK are recorded as features x1, x2, and x3, respectively.

[0163] 2. Extraction and statistics of offspring sequence feature information (species characteristics)

[0164] Using known simulated data, GC content and oligonucleotide identifiers were statistically extracted as non-family bioinformatics features to assist in sequence differentiation. For global statistical information such as GC content, read length features were sorted out, the functional relationship between species discrimination and read length was calculated, and the minimum length threshold for tolerable erroneous bases was determined. For local base information such as oligonucleotide identifiers, the species differentiation efficiency and accuracy of multiple short oligonucleotide lengths were tested according to the species genome size to reduce the sensitivity to erroneous bases.

[0165] 1. Count the GC content of read length

[0166] GC content refers to the percentage of the sum of the number of guanine and cytosine in a DNA sequence to the total number of bases in the sequence. The GC content in each long read sequence of the three simulated data was statistically calculated and recorded as feature x4. The statistical features of the GC content of the three simulated data are as follows: Figure 2 shown.

[0167] 2. Oligonucleotide identification for statistical read length

[0168] In this scheme, the frequency of oligonucleotide identification based on long read length is counted. The types of oligonucleotide identification units are 4 k When k is 3, extract each long read length 3-mer fragment unit with a length of 3bp, and count the frequency of these fragment units. The number of oligonucleotide identification units is 4 to the third power, resulting in 64 types of fragment units; in addition, due to the double helix structure of DNA and the characteristics of base complementary pairing, in actual sequencing, the two complementary chains cannot predict whether the sequencing chain is the template chain. In actual statistics, the two complementary fragment units are combined and counted, that is, a total of 32 types of fragment units are obtained. The frequency of occurrence of these fragment units in the read length is counted to obtain the statistics of 32 dimensions of oligonucleotide identification, recorded as features x5, x6, ..., x36. The 32-dimensional data feature matrix is ​​a sparse matrix.

[0169] 3. Marking high-confidence read lengths based on family feature information

[0170] Use family feature information to mark high-confidence read sequences in the offspring for statistical analysis and clustering of long read sequences in the offspring. The specific approach is to extract fragment units (K-mer units or Strobemer units) in the offspring that contain a considerable number of paternal or maternal-specific and parental-shared fragment units, that is, features x1, x2, and x3 must simultaneously exceed the corresponding thresholds, which can be marked as high-confidence read sequences in the offspring. In addition, features x1 and x2 can also be used to classify single-molecule long reads of the offspring as paternal or maternal-specific, which can lay the foundation for the next step of fully typing the haplotype genome assembly of the host species.

[0171] Since K-mer performs poorly in sequencing with high error rates, its features are buried by high error rates or random short reads, resulting in some host sequences being likely to be commensal or contaminated sequences. Further unsupervised clustering is required through other x4-x36 feature combination systems to further improve the sensitivity of host sequence recognition.

[0172] The numbers of highly confident reads marked in the progeny data of PacBio HiFi simulation data, PacBio CLR simulation data, and ONT simulation data were 285,506, 276,530, and 95,295, respectively.

[0173] 4. Multidimensional Information Clustering

[0174] In view of the different distribution characteristics of the read lengths of the three generations on three different platforms, all features were normalized. At the same time, using the K-mer marking feature, unsupervised multidimensional information clustering was performed to obtain the sequence of the offspring in the data of different platforms marked under the K-mer unit.

[0175] 1. Unsupervised Clustering

[0176] After obtaining the feature matrix of x1-x36, all long reads were clustered using the Bayesian-Gaussian mixture model. Specifically, for feature diversification, linear dimensionality reduction principal component analysis was used. When processing high-dimensional complex information space, the original feature variables were transformed to form independent principal components, objectively and truly capturing the elements with higher contribution. According to the data density distribution characteristics after dimensionality reduction, the Gaussian probability density function (normal distribution curve) was used to accurately quantify the principal components, and the Bayesian-Gaussian mixture model was applied for clustering.

[0177] The algorithm uses the initial unlabeled data to solve the prior probability of the parameters of each Gaussian model; then it uses the maximum likelihood function to iteratively solve, and uses the Bayesian method to find the maximum posterior probability class. The algorithm can effectively reduce the interference of abnormal data noise points caused by high error rates; reduce dependence on initial parameters, and benefit the development of unknown species projects.

[0178] 2. Target read length group identification based on high-confidence host sequences

[0179] Since the core of the algorithm uses an unsupervised clustering algorithm that does not introduce any prior knowledge, one of the main issues is how to mark the clusters of the clustering results. For the family genetic model, features x1-x3 should only be enriched in the host species reads, and this can be used to further identify the host species sequence clusters for all unsupervised clustering results.

[0180] Therefore, in actual operation, the density distribution of features x1-x3 corresponding to all long read sequences in the read groups of each cluster is first counted, and then the high-confidence read sequences in each read group are marked according to the actual thresholds of x1>0.1%, x2>0.1%, and x3>0.5%. All read groups are sorted according to the ratio of the number of high-confidence read sequences contained in them to the number of all read sequences in the read group, and the top two read groups (the best or second best cluster group) are taken as the target read groups for the host organism to be selected.

[0181] 3. Integration based on cyclic clustering results

[0182] Due to the randomness of unsupervised clustering algorithms and the complexity of sequence characteristics affected by sequencing errors, some host long reads cannot be guaranteed to be enriched in the same or several result groups every time, affecting the integrity and purity of the host species clustering group, and thus affecting the effect and efficiency of the final sequence decontamination.

[0183] In actual implementation, 10 cycles were used to perform the previous unsupervised clustering on the same data set, and then the target read groups for each label were integrated. After 10 cycles, 20 best or second-best cluster groups were finally obtained. The frequency of each long read sequence appearing in the 20 best or second-best cluster groups was counted, and the long read sequences with an appearance frequency of more than 80% (more than 16 times) were selected and integrated into the final host organism read group, while ensuring the integrity and reliability of the screening results.

[0184] V. Result Evaluation

[0185] Through the results of family and model-based high-dimensional clustering, each read of PacBio HiFi simulated data, PacBio CLR simulated data, and ONT simulated data can be assigned a classification label. The classification labels given in advance in the simulated data were used as the true set for evaluation to evaluate the accuracy of the decontamination of this method. At the same time, the same evaluation was performed on the mainstream decontamination tools of the same type on the market using the same set of data, and finally the accuracy, recall rate, F1-score and other indicators of the three data in different analysis and classification methods were obtained. For non-parameter tools, most non-parameter tools cannot handle long read sequences with more than 20,000 entries. Referring to Table 2, in comparison with the benchmark tool MetaBBCC-LR, the two modes of the present invention (K-mer and Strombmer) both achieved better performance in three simulated data, with each index greater than 90%, among which the Strombmer algorithm accuracy, recall rate and F1-score were all above 97%; for the reference tools Centrifuge, Kraken2 and MetaMaps that rely on large databases, both PacBio HiFi simulation data and PacBio CLR simulation data were close to 100%. However, in the ONT simulation data with a higher error rate, the performance declined. Based on the above results, it can be seen that the two modes of the present invention do not rely on any database and can obtain similar high-quality decontamination effects.

[0186] Table 2. Comparison of the performance of different sequence classification software and the method of the present invention on simulated data

[0187]

[0188] 6. Genome assembly using decontaminated host data

[0189] 1. Host genome assembly

[0190] For the long read data of the purified host, Canu or Flye are used to perform high-quality genome assembly. With the reference genome sequence as the evaluation criteria, QUAST is used to comprehensively and objectively evaluate the integrity, length, continuity and other aspects of the genome assembly sequence before and after decontamination. As a result, in this embodiment, in all PacBioHiFi, PacBio CLR and ONT data, the total length of assembly (Total length) of Canu or Flye any assembly software after implementing the decontamination process is greatly reduced to the reference genome length range (Reference length); At the same time, the sequence length (Unaligned length) not aligned to the reference genome sequence is significantly reduced. This shows that after the three generations of long reads of decontamination (symbiosis), the remaining pure long reads cannot support the genome assembly of other mixed species, effectively eliminating the interference of other alien species. The high recall rate of the host long read ensures the assembly quality of the host genome itself. In addition, due to the reduction of the impact of the pollution (symbiosis) sequence, the genome assembly map is greatly simplified, which greatly reduces the difficulty and complexity of the genome assembly process based on graph theory. This is reflected in the fact that the indicators such as genome assembly errors (misassemblies) and assembly mismatches (mismatch) of most of the data in the table are reduced. Especially for the PacBio HiFi simulation data with the highest single-base sequencing accuracy, the improvement effect is most obvious. Affected by this, the continuity (N50) of the genome assembly of most of the data is also improved. For example, in the ONT simulation data with the longest sequencing read length, although its single-base sequencing error rate is the highest, both assembly softwares increase the assembly sequence N50 to about 13M. This fully illustrates that the contaminated sequence filtering method described in the present invention can efficiently and stably obtain pure host genomes, and is not affected by sequencing platforms and assembly software.

[0191] Table 3. Referenced evaluation of host genome assembly quality before and after removing contaminating sequences using the present invention

[0192]

[0193] 2. Completeness and contamination of host genome assembly

[0194] Compare the length, GC content, integrity, contamination and other indicators of the host genome assembly before and after the contamination sequence is filtered. The results refer to Table 4, which further illustrates the improvement of the quality of the host genome assembly by decontamination. In this embodiment, in the simulated sequencing data of all PacBio HiFi, PacBio CLR and ONT, using any assembly software of Canu or Flye, after the decontamination process is implemented, the assembly length and GC content are closer to the reference genome sequence, which is consistent with the QUAST evaluation. In addition, each group of data has greatly reduced the contamination on the basis of keeping the assembly integrity almost unchanged, from more than 350% to less than 2%. This fully illustrates that the contamination sequence filtering method described in the present invention can effectively remove the genome sequence of foreign species in the mixed data, reduce the genome contamination, and is not affected by the sequencing platform and assembly software. In addition, it is worth mentioning that because the accuracy and recall rate of decontamination are maintained at a high level, the genome assembly GC content obtained by removing the contamination sequence using the present invention is closer to the actual content of the reference genome of the species itself.

[0195] Table 4. Host genome assembly integrity and contamination evaluation before and after removing contamination sequences from three types of simulated data using the present invention

[0196]

[0197] Example 2: Real data of long-read sequencing of three platforms of host-symbiont based on family

[0198] In order to better illustrate the practicality and innovation of the present invention, real samples with family information were used to construct DNA long-read sequencing data with host-contaminant-symbiont as the sample source, reflecting the performance of the present invention in practical applications.

[0199] In this example, three samples with the paternal number HG003, the maternal number HG004, and the offspring number HG002 in the 1000 Genomes Database were also used as the parental source families of the family samples. The NGS data of the paternal parent extracted from the GIAB public database was 44,485,998 paired-end (PE) sequences, and the NGS data of the maternal parent was 46,384,012 paired-end sequences.

[0200] Real PacBio RSIICLR platform data: human chromosome 19 of HG002 as the host, two bacteria (Bacillus subtilis and Lactobacillus fermentum) from the microbial dataset (Mock10) of the ZymoBIOOMICS microbiome standard as symbionts with the host, and eight other bacteria and yeast as random contaminants. In this case, chimpanzee chromosome 21 was selected as a symbiont with highly similar sequences to test the sensitivity to sequence similarity between host and symbiont species. Human HG002 raw sequencing data was downloaded from NCBI, and raw sequencing data of chromosome 19 with 50× alignment to the HG002 reference genome (GCA_011064465.1) was extracted (Shumate et al., 2020). Chimpanzee raw sequencing data was downloaded (Logsdon et al., 2021), and sequencing data with 50× alignment to the chimpanzee chromosome 21 reference genome (GCA_000001515.5) was extracted. The microbial raw data of Mock10 were downloaded from the NCBI website (McIntyre et al., 2019). The same ratio of 10 bacteria and yeast was maintained in the mixed data.

[0201] Real ONT PromethION platform data: The composition is the same as the real PacBio RSIICLR data. Human HG002 raw sequencing data was downloaded from NCBI, and raw sequencing data of chromosome 19 was extracted with 50× alignment to the HG002 reference genome (GCA_011064465.1) (Shumate et al., 2020). Chimpanzee raw sequencing data was downloaded (Logsdon et al., 2021), and sequencing data of chromosome 21 of chimpanzee was extracted with 50× alignment to the reference genome (GCA_000001515.5). The same ratio of 10 bacteria and yeast was maintained in the mixed data. Mock10 microbial raw data was downloaded from the NCBI website (Nicholls et al., 2019).

[0202] Real Cyclone WT platform data: The long read data of HG002 was obtained using the long read sequencing technology platform of the Cyclone WT platform, a domestic sequencing platform of Shenzhen BGI Life Sciences Research Institute, and a total of 836,557 reads were extracted and aligned to chromosome 19 of the HG002 reference genome (GCA_011064465.1). The Mock10 data with a total of 1,420,920 reads were obtained using ZymoBIOOMICS microbial community standard sequencing. Two bacteria (Bacillus subtilis and Lactobacillus fermentum) were used as symbiotic bacteria with the host, and the other eight bacteria and yeast were used as random contaminants. In this embodiment, not only the most mainstream single-molecule sequencing technology platforms on the market are covered, but also the newly developed new generation single-molecule platform Cyclone WT with complete domestic intellectual property rights. Detailed information on the artificially mixed species components of the real data of the above three long-read sequencing platforms is shown in Table 5.

[0203] Table 5. PacBio RSIICLR, ONT PromethION and Cyclone WT platform real sequencing data composition information table

[0204]

[0205] 1. Extraction of family genetic characteristic information (characteristic fragment set)

[0206] 1. Extraction of family genetic feature sequence information based on K-mer units

[0207] When constructing the genetic characteristics information of the family, the K-mer generation tool (meryl tool v1.3) was first used to perform K-mer characterization on the NGS data of the paternal HG003 and maternal HG004, quickly generate K-mer information, and perform set operations to generate the paternal-specific set POK, the maternal-specific set MOK, and the parent-to-parent common set SK. The specific process is as follows: the meryl tool is used to extract the 21-mer fragment unit of each read length in the paternal and maternal sequencing data, that is, firstly, based on the sequence fragment unit with a fixed length of 21 bases, each long read length of the paternal and maternal is extracted from the first base in the form of a sliding window (step size of 1), for example, to sequentially obtain the 21-mer fragment units of all read lengths. Secondly, the 21-mer sequence frequencies of the paternal and maternal are counted respectively, and the sequencing error peaks and repeated sequence peaks are deleted using the frequency statistical distribution, that is, the ultra-low frequency and ultra-high frequency 21-mer sequences are deleted. At this time, the paternal 21-mer set is recorded as set P, and the maternal 21-mer set is recorded as set M. Then, in order to obtain the paternal-specific 21-mer set, set P is used to take the difference set of set M, and finally the paternal-specific set POK is obtained. In order to obtain the maternal-specific 21-mer set, set M is used to take the difference set of set P, and finally the maternal-specific set MOK is obtained. Finally, the intersection of set P and set M is taken, which is the paternal and maternal common 21-mer set, recorded as SK.

[0208] Finally, the characteristic sequence information was obtained to establish a family characteristic set with a K-mer fragment unit of 21bp in length, with a total of 14,042,462 21-mers in the paternal-specific set POK, 16,747,856 21-mers in the maternal-specific set MOK, and 76,947,365 21-mers in the parent-common set SK. The total number of 21-mers with family characteristics was 107,737,683.

[0209] 2. Extraction of characteristic information of pedigree genetic sequences based on strobemer unit fragments

[0210] Use Strobemer fragment units to extract family characteristic sequence information. This scheme chooses to use random Strobes with n=2, w_min=15, and w_max=30 as a practical technology. First, obtain a 40-mer K-mer feature database. Secondly, convert the 40-mers one by one into corresponding Strobemer representative sequences. At this time, the paternal Strobemer set is recorded as set P, and the maternal Strobemer set is recorded as set M. Then, in order to obtain the paternal-specific Strobemer set, use set P to take the difference set of set M, and finally obtain the paternal-specific set POK. In order to obtain the maternal-specific Strobemer set, use set M to take the difference set of set P, and finally obtain the maternal-specific set MOK. Finally, the intersection of set P and set M is the paternal and maternal common Strobemer set, recorded as SK. Finally, the parent-specific Strobemer sets POK, MOK and the common Strobemer set SK are obtained.

[0211] Finally, the characteristic sequence information obtained was a family characteristic set of Strobemer fragment units with a length of 40bp, with a total of 35,106,816 Strobemers in the paternal-specific set POK, 39,174,692 Strobemers in the maternal-specific set MOK, and 148,584,406 Strobemers in the parent-common set SK. The total number of family characteristic Strobemers was 222,865,914.

[0212] In the extraction of family characteristics, the paternal-specific set POK, the maternal-specific set MOK, and the parent-common set SK are recorded as features x1, x2, and x3, respectively.

[0213] 2. Extraction and statistics of offspring sequence feature information (species characteristics)

[0214] 1. Count the GC content of read length

[0215] The GC content in the long read sequences of the PacBio RSIICLR dataset, the ONT PromethION platform dataset, and the CycloneWT platform dataset was statistically calculated and recorded as feature x4. The statistical features of the GC content of the three real sequencing data are as follows: Figure 3 shown.

[0216] 2. Oligonucleotide identification for statistical read length

[0217] In this scheme, the frequency of oligonucleotide identification based on long read length is counted. The types of oligonucleotide identification units are 4 kWhen k is 3, extract each long read length 3-mer fragment unit with a length of 3bp, and count the frequency of these fragment units. The number of oligonucleotide identification units is 4 to the third power, resulting in 64 types of fragment units; in addition, due to the double helix structure of DNA and the characteristics of base complementary pairing, in actual sequencing, it is impossible to predict whether the two complementary chains are template chains. In actual statistics, the two complementary fragment units are combined and counted, that is, a total of 32 types of fragment units are obtained. The frequency of occurrence of these fragment units in the read length is counted to obtain the statistics of 32 dimensions of oligonucleotide identification, recorded as features x5, x6, ..., x36. The 32-dimensional data belongs to a sparse matrix.

[0218] 3. Marking high-confidence read lengths based on family feature information

[0219] The high-confidence read sequences in the offspring are marked with family feature information for statistical analysis and clustering of the offspring long read sequences of the PacBio RSIICLR dataset, ONT PromethION platform dataset, and Cyclone WT platform dataset. The specific approach is to extract the offspring that contain a considerable number of fragment units (K-mer units or Strobemer units) that are unique to the father or mother and common to the parents, that is, the features x1, x2, and x3 must exceed the corresponding thresholds at the same time, and they can be marked as high-confidence read sequences in the offspring. The thresholds set for x1, x2, and x3 are 0.1%, 0.1%, and 0.5%, respectively. In addition, the features x1 and x2 can also be used to classify the single-molecule long reads of the offspring as paternal or maternal specific, which can provide a good data basis for the next step of fully typing the haplotype genome assembly of the host species.

[0220] The remaining sequences may contain host sequences buried by commensal organisms, contaminated sequences, high error rates, or random short reads, and further unsupervised clustering is required through other x4-x36 feature combination systems to further improve the sensitivity of host sequence identification.

[0221] 4. Multidimensional Information Clustering

[0222] In view of the different distribution characteristics of the read lengths of the long read data of several platforms, all features are normalized to the same length. At the same time, the K-mer marker feature is used to further integrate the application of biological information and experimental technology features in the pollution filtering system.

[0223] 1. Unsupervised Clustering

[0224] After obtaining the feature matrix of x1-x36, all long reads were clustered using the Bayesian-Gaussian mixture model. Specifically, for feature diversification, linear dimensionality reduction principal component analysis was used to transform the original feature variables to form independent principal components. According to the data density distribution characteristics after dimensionality reduction, the Gaussian probability density function (normal distribution curve) was used to accurately quantify the principal components, and the Bayesian-Gaussian mixture model was applied for clustering.

[0225] 2. Target class identification based on high-confidence host sequences

[0226] First, the density distribution of features x1-x3 corresponding to all long read sequences in the read groups of each cluster is counted, and then the high-confidence read sequences in each read group are marked according to the actual thresholds of x1>0.1%, x2>0.1%, and x3>0.5%. All read groups are sorted according to the ratio of the number of high-confidence read sequences contained in them to the number of all read sequences in the read group, and the top two read groups (the best or second best cluster group) are taken as the target read groups for the host organism to be selected.

[0227] 3. Integration based on cyclic clustering results

[0228] The same data set was clustered unsupervisedly for 10 cycles, and the target read groups for each marker were integrated. After 10 cycles, 20 best or second-best clusters were finally obtained. The frequency of each long read sequence appearing in the 20 best or second-best clusters was counted, and the long read sequences with an appearance frequency of more than 80% (more than 16 times) were selected and integrated into the final host organism read group, while ensuring the integrity and reliability of the screening results.

[0229] The number of hosts finally marked in the PacBio RSIICLR platform dataset, ONT PromethION platform dataset, and CycloneWT platform dataset were 358,368, 325,894, and 166,695, respectively.

[0230] V. Result Evaluation

[0231] Through the results of family and model-based high-dimensional clustering, each read length of the PacBio RSIICLR platform data set, ONTPromethION platform data set, and Cyclone WT platform data set was assigned a classification label, and the accuracy and efficiency of the decontamination of this method were evaluated based on the classification labels. At the same time, the effects of similar tools in read length classification were also performed here, including MetaBCC-LR, Centrifuge, Kraken2, and MetaMaps. It should be pointed out that MetaBCC-LR is a classification tool that does not require a reference database, while the others are reference classification tools based on existing reference genome libraries. Most of the other commonly used non-reference read length classification tools on the market cannot handle long reads of more than 20,000, so they are not shown here. Referring to Table 6, in the data performance of the PacBio RSIICLR platform, the accuracy of the two marker family feature algorithms, Strobemer and K-mer, proposed in the present invention are 99.3% and 98.2%, respectively, which are higher than all other tools; in terms of the recall rate index, they are generally at a high level, 94.1% and 96.2%, respectively, which are much higher than the 56.70% of the non-reference classification tool Meta BCC-LR, and slightly lower than the reference classification tools Centrifuge and Kraken2; in addition, in terms of the F1-score index, the present invention also obtains the top ranking performance, which are 96.6% and 97.2%, respectively. In the data performance of the ONT PromethION platform, the specific performance of the two marker family feature algorithms, Strobemer and K-mer, proposed in the present invention is that the Strobemer algorithm is equivalent to Centrifuge and Kraken2; whether it is the recall rate or the F1-score, the results of this method are better than the similar non-reference tool MetaBCC-LR. In the real sequencing data of ONT, the reference tool is generally slightly better than the non-reference tool. In terms of data performance on the Cyclone WT platform, the two marker family feature algorithms, Strobemer and K-mer, proposed by the present invention, have accuracy rates of 98.78% and 98.69%, respectively, which are good results; in terms of recall rate indicators, they are generally at a high level, 99.85% and 97.71%, respectively, which are much higher than similar non-parameter tools MetaBCC-LR; only slightly lower than the parameter classification tools Centrifuge and Kraken2. In addition, in terms of F1-score indicators, the present invention also achieved the top ranking performance, 99.31% and 98.20%, respectively, which is nearly twice the performance of the non-parameter tool MetaBCC-LR.

[0232] In general, the algorithm model for removing contaminated sequences proposed in the present invention has good performance on data from three different long-read platforms in two marker family feature modes (Strobemer and K-mer). Overall, various indicators are stable, and most indicators are better than the non-reference tool MetaBCC-LR software, among which the F1-score is nearly doubled. It is also worth mentioning that in the two family feature information extraction modes proposed in the present invention, Strobemer performs better than the K-mer algorithm, both reaching more than 91%. This is comparable to the results of MetaMaps, the best performing tool with references, and is more stable than the widely used reference tool Kraken2. Therefore, the use of the Strobemer algorithm in the present invention can obtain the purest data with the best performance, and perform downstream assembly analysis based on this.

[0233] Table 6. Comparison of the performance of different sequence classification software and the method of the present invention on real data of three platforms

[0234]

[0235] 6. Genome assembly analysis using decontaminated host data

[0236] 1. Host genome assembly and reference evaluation

[0237] For the long read length data of the purified host, Canu or Flye is used to perform high-quality genome assembly. In this embodiment, the reference genome sequence is used as the evaluation criterion, and QUAST is used to evaluate the integrity, length, continuity and other aspects of the genome assembly sequence before and after decontamination. The results refer to Table 7. After the decontamination process is implemented, the data of the three platforms are closer to the reference genome length (Referencelength) for the total length of the host genome assembly; At the same time, the sequence length (Unaligned length) that is not aligned to the reference genome sequence also decreases significantly. In addition, after excluding the interference of the alien species sequence on the genome assembly map, the negative indicators such as genome assembly errors (misassemblies) and assembly mismatches (mismatch) of most of the data have decreased. At the same time, the continuity (N50) of the genome assembly has been improved. This shows that for the real sequencing data of different platforms, the contaminated sequence filtering method described in this patent can still efficiently and stably obtain the pure host genome assembly result.

[0238] Table 7. Referenced evaluation of the host genome assembly quality before and after removing contamination sequences from real data of three platforms by the present invention

[0239]

[0240] 2. Host genome assembly completeness and contamination

[0241] In order to more comprehensively and quantitatively analyze the impact of decontamination on host genome assembly, this embodiment also uses genome-specific K-mers of various species (including host, symbiotic species, and random contamination species) (i.e., K-mers that exist exclusively in the genome of this species and not in the genome of any other species) to evaluate the integrity and contamination of genome assembly results based on real sequencing data of PacBio RSIICLR platform, ONTPromethION platform, and Cyclone WT platform. The results refer to Table 8. The integrity of the genome assembly based on the filtered host genome data is maintained well, while the contamination decreases by 45.7 times, 11.9 times, and 20.7 times, respectively. The contamination is significantly improved, which will provide more accurate genome data guarantee for the analysis of downstream applications.

[0242] Table 8. Evaluation of the completeness and contamination of the host genome assembly before and after removing contamination sequences from three real data using the present invention

[0243]

[0244] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be equivalent replacement methods and are included in the protection scope of the present invention.

Claims

1. A method for filtering sequencing data of an organism, characterized in that: The following steps are involved: S110: clustering the sequencing data of the organism according to the content and species characteristics of the characteristic fragment set to obtain a plurality of read length groups, wherein the characteristic fragment set includes at least one of a parental common characteristic fragment set, a paternal unique characteristic fragment set, and a maternal unique characteristic fragment set; S120: marking a target read length group having a high-confidence read length sequence from the plurality of read length groups, wherein the high-confidence read length sequence is a sequence whose content of the characteristic fragment set is above a threshold value; S130: Repeat the cycle of S110 and S120 to screen out read length sequences whose occurrence frequency in all target read length groups is higher than a set value.

2. The filtering method according to claim 1, characterized in that: The characteristic fragment set is obtained by performing set operation based on the paternal Kmer sequence set and the maternal Kmer sequence set; Preferably, the set operation includes taking intersection and difference; Preferably, the characteristic fragment set obtained by performing set operation based on the paternal Kmer sequence set and the maternal Kmer sequence set includes: Based on the sequencing data of the paternal and maternal copies of the organism, a paternal Kmer sequence set and a maternal Kmer sequence set are obtained; Performing a set operation on the paternal Kmer sequence set and the maternal Kmer sequence set to obtain a characteristic fragment set; Preferably, the characteristic fragment set is a Kmer sequence set; Preferably, the length of the Kmer sequence in the Kmer sequence set is 10 to 100 bases; Preferably, before the set operation, the method further includes using the Strobemer method to convert the paternal Kmer sequence set and the maternal Kmer sequence set into a paternal Strobemer sequence set and a maternal Strobemer sequence set; the characteristic fragment set is a Strobemer sequence set; Preferably, the Strobemer method includes at least one of minstrobes, hybridstrobes, randstrobes, mixedstrobes, altstrobes, and multistrobes; Preferably, the strobemer sequence in the strobemer sequence set is based on a sliding window of variable size, each window extracts an l-mer, and connects n consecutive l-mers head to tail to obtain n*l-mer; wherein l is 10 to 100, n is 2 to 10, and the variable size has a minimum value w min and the maximum value w max , w min 5~20,w max is 10 to 50, and w min <w max ; Preferably, the sequencing data of the paternal and maternal copies of the organism are second generation sequencing data; Preferably, the sequencing data of the organism is third-generation sequencing data; Preferably, the third-generation sequencing data is any one of PacBio HiFi, PacBio CLR, ONT, and Cyclone WT sequencing data.

3. The filtering method according to claim 1, characterized in that: In S120, the sequence whose content of the characteristic segment set is above the threshold is a sequence whose content of the paternal-specific characteristic segment set is above the first threshold, whose content of the maternal-specific characteristic segment set is above the second threshold, and whose content of the parent-common characteristic segment set is above the third threshold; Preferably, the first threshold, the second threshold and the third threshold are independently 0.01 to 1%; Preferably, the first threshold is independently 0.05-0.5%, the second threshold is independently 0.05-0.5%, and the third threshold is independently 0.2-0.8%; Preferably, the first threshold is independently 0.05-0.2%, the second threshold is independently 0.05-0.2%, and the third threshold is independently 0.4-0.6%.

4. The filtering method according to claim 1, characterized in that: The species characteristics include at least one of GC content and frequency of oligonucleotide markers; Preferably, the length of the oligonucleotide marker is 2 to 10; Preferably, the length of the oligonucleotide marker is 2 to 4; Preferably, the length of the oligonucleotide marker is 3; Preferably, among the frequencies of the oligonucleotide markers, the frequencies of the complementary oligonucleotides are combined for calculation.

5. The filtering method according to claim 1, characterized in that: In S130, the number of cycles is 5 to 20 times, and / or the set value is 60 to 99%.

6. A method for assembling sequencing data of an organism, characterized in that: The following steps are involved: Filtering according to the filtering method according to any one of claims 1 to 5, and then assembling the filtered read sequences; Preferably, the assembly tool is selected from Canu and Flye; Preferably, it also includes classifying the read sequences into paternal-specific read sequences and maternal-specific read sequences according to the paternal-specific characteristic fragment set and the maternal-specific characteristic fragment set, and assembling them separately.

7. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer-executable instructions, and the computer-executable instructions are used to enable a computer to execute the filtering method described in any one of claims 1 to 5 or the assembling method described in claim 6.

8. An electronic device, characterized in that: The electronic device includes a processor and a memory, wherein the memory stores a computer program that can be run on the processor, and the processor implements the filtering method described in any one of claims 1 to 5 or the assembly method described in claim 6 when running the computer program.

9. A filtering system for sequencing data, characterized in that: include: A clustering module, wherein the clustering module is used to cluster the sequencing data of the organism according to the content and species characteristics of the characteristic fragment set to obtain a plurality of read length groups, wherein the characteristic fragment set includes at least one of a parental common characteristic fragment set, a paternal unique characteristic fragment set, and a maternal unique characteristic fragment set; A marking module, the marking module is used to mark a target read length group having a high-confidence read length sequence from the plurality of read length groups, the high-confidence read length sequence being a sequence whose content of the characteristic fragment set is above a threshold value; A screening module is used to screen out read sequences whose occurrence frequency in all target read length groups is higher than a set value after the clustering module and the marking module are repeatedly cycled.

10. A system for assembling sequencing data, characterized in that: include: The filtration system according to claim 9; An assembly module is used to assemble the read length sequences screened by the filtering system.

Images