A phage metagenomic sequencing data host sequence filtering system and method

Abstract

The application belongs to the technical field of data filtering, and discloses a bacteriophage metagenomic sequencing data host sequence filtering method; comprising: obtaining bacteriophage metagenomic sequencing data of a sample to be analyzed and performing pretreatment to obtain quality filtered read data; calculating a host abundance background value corresponding to the sample to be analyzed, configuring a corresponding host genome library, and setting an alignment determination threshold; performing source classification, outputting a host mapping read set, a bacteriophage candidate read set, and a boundary candidate read set; performing joint source re-determination to identify remedial fragments belonging to non-host sequences; performing sequence assembly to obtain candidate contigs and performing true source determination, and performing reverse correction based on the determination result; identifying a deviation state and updating and compensating the alignment determination threshold of the sample to be analyzed in a subsequent batch; and improving the accuracy of bacteriophage metagenomic sequencing data host sequence filtering.

Description

Technical Field

[0001] This application relates to the field of data filtering technology, and more specifically, to a host sequence filtering system and method for phage metagenomic sequencing data. Background Technology

[0002] Accurate analysis of phage metagenomic sequencing data is a crucial step in research fields such as environmental microbiome and virology, directly impacting the understanding of phage diversity and its interaction mechanisms with the host. Existing traditional filtering methods typically use alignment software to map raw sequencing data to the host reference genome. While these methods can perform basic sequence cleaning under ideal conditions with simple host backgrounds and short sequencing data, they fail to achieve the desired results in complex environmental interference scenarios and still face numerous technical bottlenecks.

[0003] In real-world metagenomic sequencing scenarios, phage sequences not only exist in a free state but may also integrate into the host genome as prophages. However, traditional filtering methods lack effective segmentation of the host-phage boundary region, easily deleting entire reads spanning host and non-host regions as host sequences, leading to the loss of phage integration boundary information. Furthermore, when processing paired-end sequencing data, traditional filtering methods can only determine the alignment status of individual reads, failing to combine insert length and orientation information to extract pairing conflict features and perform joint origin determination. This easily leads to the fragmentation of structural variation features with pairing relationships, resulting in the misjudgment and removal of valuable non-host sequences. Simultaneously, alignments at the short read level have a very high risk of missed detection, and traditional filtering methods tend to... Traditional methods assume that filtered reads can be directly used for subsequent analysis. This leads to the generation of erroneous sequences in the assembly stage due to missed host sequences in the initial stage. The lack of an error correction mechanism to trace back from long sequence contigs to the underlying reads makes it impossible to identify and remove contaminated fragments in the final assembly. On the other hand, different types of samples have different host abundance backgrounds. Traditional filtering methods usually use fixed alignment thresholds and do not consider adaptive feedback based on read migration deviations. This results in a large number of contaminated fragments escaping when the host background intensity is high, or excessive removal of low-abundance phages when the host background intensity is low. This seriously affects the generalization ability of the filtering system in variable and complex sample environments and the reliability of the final sequencing analysis results.

[0004] In view of this, this application proposes a host sequence filtering system and method for phage metagenomic sequencing data to solve the above problems. Summary of the Invention

[0005] To overcome the aforementioned deficiencies of the prior art and to achieve the above objectives, this application provides the following technical solution: a method for filtering host sequences in phage metagenomic sequencing data, comprising:

[0006] S1. Obtain the phage metagenomic sequencing data of the sample to be analyzed and preprocess it to obtain quality filtered read data; extract the compositional features of the quality filtered read data and calculate the corresponding host abundance background value of the sample to be analyzed, while configuring the corresponding host genome library and setting the comparison judgment threshold.

[0007] S2. Map and align the quality-filtered read data to the host genome database, extract the read mapping status features of each read in the mapping and alignment results, perform source classification on the quality-filtered read data in combination with the alignment judgment threshold, and output the host mapped read set, the phage candidate read set, and the boundary candidate read set;

[0008] S3. Extract local breakpoints and pairing conflict features from the boundary candidate read set, segment the fragment sequence based on the local breakpoints, perform joint source re-determination on the segmentation results in combination with the pairing conflict features, identify the rescue fragment that belongs to the non-host sequence, and transfer the rescue fragment into the phage candidate read set;

[0009] S4. Assemble the sequence of the candidate reads of the transferred phage to obtain candidate contigs; determine the true source of the candidate contigs, establish a mapping link between the candidate contigs and the corresponding supporting reads, and make reverse corrections to the source classification of the corresponding supporting reads based on the determination results.

[0010] S5. Based on the reverse correction results, separate host-origin reads from the quality-filtered read data and output the filtered phage metagenomic sequencing data; calculate the read migration rate during the reverse correction process, and update and compensate the comparison threshold for subsequent batches of samples to be analyzed based on the deviation between the read migration rate and the corresponding host abundance background value.

[0011] Furthermore, methods for extracting compositional features from quality-filtered read data and calculating the corresponding host abundance background value for the sample to be analyzed include:

[0012] Sequence sampling is performed on the quality filtered read data to identify the percentage of specific bases and fixed-length nucleotide short strings in each sampled read, which are then combined to form compositional features.

[0013] A fixed-length nucleotide short string is input into a preset host feature hash index for matching, and the number of read segments that match is counted.

[0014] The initial ratio is obtained by calculating the ratio of the number of matched reads to the total number of corresponding sampled reads; the overall base content distribution curve is plotted based on the specific base content percentage of each sampled read, the overlap area between the overall base content distribution curve and the corresponding base distribution interval is identified, and the base overlap is output; the initial ratio is weighted and corrected using the base overlap to obtain the host abundance background value of the corresponding sample to be analyzed.

[0015] Furthermore, the methods for extracting the reading segment mapping status features of each reading segment in the mapping comparison results and combining them with the comparison judgment threshold to perform source classification on the quality-filtered reading segment data include:

[0016] The genome sequence in the host genome library is identified as a reference index. The quality filtered read data is matched with the reference index to generate sequence alignment records as mapping comparison results.

[0017] The sequence coverage ratio, matching similarity, sequence alignment features with soft shearing markers, and position features of the paired read segments in the parsing mapping and alignment results are combined into read segment mapping state features.

[0018] Reads with sequence coverage and matching similarity not lower than the corresponding alignment threshold are added to the host-mapped read set; reads that cannot be mapped to the host genome library at either end are added to the phage candidate read set.

[0019] Unmapped reads with soft shearing markers in the sequence alignment features, and reads with one end mapped and the other end unmapped or of abnormal length in the position features of the paired reads, are added to the boundary candidate read set.

[0020] Furthermore, methods for segmenting sequences based on local breakpoints include:

[0021] Traverse the sequence alignment records corresponding to each read in the boundary candidate read set and identify the soft cut markers; at the same time, extract the alignment direction information and insertion fragment length information of the paired reads and combine them as pairing conflict features.

[0022] Identify the cleaved base sequence indicated by the soft shearing marker and determine whether the number of bases in the cleaved base sequence reaches the preset truncation length lower limit; if the number of bases reaches the preset truncation length lower limit, then confirm the specific boundary coordinates of the corresponding read segment as a local breakpoint based on the sequence alignment record;

[0023] Based on local breakpoints, the main body of the corresponding read segment is physically truncated, the unmapped short sequences are separated, and they are output as independent segments.

[0024] Furthermore, the methods for performing joint source re-determination on the segmentation results include:

[0025] Identify the paired end read segments corresponding to the independent segmentation segments, and obtain the paired end alignment direction information and paired end inserted segment length information from the paired end conflict features corresponding to the paired end read segments;

[0026] When the paired end read is not mapped to the host genome library, and the paired end alignment direction information is reversed or the length information of the paired end inserted fragment exceeds the preset normal length range, the corresponding independent segment and the paired end read will be combined into a joint target.

[0027] The sequence composition of the joint target is extracted, and the corresponding region of the sequence composition is compared with the pre-defined conserved region of phage features. When the joint target contains the pre-defined conserved region of phage features, the joint target is identified as a rescue fragment of a non-host sequence, and the rescue fragment is transferred as a whole to the phage candidate read set.

[0028] Furthermore, methods for determining the true source of candidate contiguous groups include:

[0029] The transferred salvage fragments and the original set of phage candidate reads are input into a pre-defined metagenomic assembly structure and spliced ​​to generate a candidate contiguous group with a continuous sequence structure.

[0030] Gene coding sequences extracted from candidate contigs are compared with host genome databases for long-scale homology. The proportion of host homologous sequences corresponding to the gene coding sequence is calculated. At the same time, the sequence complexity and viral motif density of unmapped regions in the corresponding candidate contigs are obtained.

[0031] When the proportion of host homologous sequences is higher than the preset contamination judgment threshold, and the sequence complexity of the unmapped region is lower than the expected complexity threshold, the corresponding candidate contiguous group will be judged as the host contamination contiguous group.

[0032] When the proportion of host homologous sequences is lower than the safety threshold and the density of viral characteristic motifs is higher than the significant density, the corresponding candidate contigs are determined to be true phage contigs.

[0033] Furthermore, the methods for establishing a mapping link between candidate contiguous groups and corresponding supporting reads, and for making reverse corrections to the source classification of the corresponding supporting reads based on the judgment results, include:

[0034] Extract the sequence identifier set of all initial short reads used to construct each candidate contiguous group, and use each initial short read as a supporting read;

[0035] Based on the affiliation relationship between the sequence identifier set and the corresponding candidate overlapping group, a mapping link is constructed from the candidate overlapping group to each corresponding supporting read segment in the corresponding sequence identifier set;

[0036] When the corresponding candidate contiguous group is determined to be a host-contaminated contiguous group, a reverse tracing operation is performed on the corresponding mapping link to lock all abnormal support reads that constitute the host-contaminated contiguous group.

[0037] The source classification labels of all abnormal support reads were uniformly changed to host classification, and all abnormal support reads were separated from the phage candidate read set after being transferred to the salvage fragment and merged and transferred to the host-mapped read set.

[0038] Furthermore, the read migration rate during the reverse correction process is statistically analyzed, and the deviation between the read migration rate and the corresponding host abundance background value is determined by the following methods:

[0039] The total number of anomalous supporting reads that are separated from the phage candidate read set after being transferred to the rescue fragment and transferred to the host-mapped read set during the statistical reverse correction process;

[0040] The read mobility is obtained by calculating the ratio of the total number of reads in the sequence to the total number of reads in the quality-filtered read data.

[0041] Extract the host abundance background value of the corresponding sample to be analyzed, and obtain the preset escape tolerance coefficient. Calculate the product of the preset escape tolerance coefficient and the host abundance background value to obtain the maximum allowable tolerance value of the corresponding sample to be analyzed. If the read migration rate is higher than the maximum allowable tolerance value, the corresponding deviation state is determined to be a positive deviation state; otherwise, the deviation state is determined to be a stable deviation state.

[0042] Furthermore, the methods for updating and compensating the comparison thresholds for subsequent batches of samples to be analyzed include:

[0043] If the deviation is positive, obtain the sample source label of the corresponding sample to be analyzed, and calculate the ratio of the difference between the corresponding read segment mobility and the maximum allowable tolerance value.

[0044] Calculate the product of the difference ratio and the preset adjustment step size to obtain the compensation value; superimpose the corresponding comparison judgment threshold and the compensation value to obtain the compensated threshold; bind the compensated threshold with the corresponding sample source label and save it; when subsequent batches with the same sample source label appear, apply the corresponding compensated threshold as the comparison judgment threshold.

[0045] A host sequence filtering system for phage metagenomic sequencing data, used to implement a method for filtering host sequences in phage metagenomic sequencing data, comprising:

[0046] The preliminary filtering module acquires and preprocesses the phage metagenomic sequencing data of the sample to be analyzed to obtain quality-filtered read data; it extracts the compositional features of the quality-filtered read data and calculates the corresponding host abundance background value of the sample to be analyzed, while configuring the corresponding host genome library and setting the alignment judgment threshold.

[0047] The source classification module maps and aligns the quality-filtered read data to the host genome database, extracts the read mapping status features of each read in the mapping and alignment results, performs source classification on the quality-filtered read data in combination with the alignment judgment threshold, and outputs the host mapped read set, the phage candidate read set, and the boundary candidate read set.

[0048] The fragment transfer module extracts local breakpoints and pairing conflict features from the boundary candidate read set, performs fragment sequence segmentation based on local breakpoints, and performs joint source re-determination on the segmentation results in combination with pairing conflict features to identify rescue fragments belonging to non-host sequences, and transfers the rescue fragments into the phage candidate read set;

[0049] The reverse correction module assembles the sequence of the transferred phage candidate read set to obtain candidate contigs; it determines the true source of the candidate contigs, establishes a mapping link between the candidate contigs and the corresponding supporting reads, and performs reverse correction on the source classification of the corresponding supporting reads based on the determination results.

[0050] The compensation update module separates host-origin reads from the quality-filtered read data based on the reverse correction results and outputs the filtered phage metagenomic sequencing data; it calculates the read migration rate during the reverse correction process and updates and compensates the comparison threshold for subsequent batches of samples to be analyzed based on the deviation of the read migration rate from the corresponding host abundance background value; the modules are connected to each other via wired and / or wireless means.

[0051] The technical effects and advantages of the host sequence filtering system and method for phage metagenomic sequencing data proposed in this application are as follows:

[0052] This paper presents a host sequence filtering method for phage metagenomic sequencing data by acquiring phage metagenomic sequencing data of samples to be analyzed and performing host sequence filtering. Compared with existing technologies, this method improves the mismatch problem caused by applying fixed filtering conditions to samples with different host backgrounds by specifically configuring the host genome library and alignment threshold. By performing fragment sequence segmentation, it avoids the overall deletion of boundary reads containing both host and non-host regions, preventing the loss of effective phage fragments. By performing joint source re-determination, it solves the problem of broken evidence chains of structural anomalies caused by splitting paired reads for independent processing. By performing true source determination and reverse correction, it enables the re-identification and recovery of host contamination sequences missed in the initial screening stage at the long sequence scale. By updating and compensating the alignment threshold of subsequent batches based on the deviation between read migration rate and host abundance background value, it improves the adaptability of the host sequence filtering process in cross-batch applications of the same type. This improves the accuracy of host sequence filtering for phage metagenomic sequencing data and enhances the stability and generalization ability of the filtering results under different sample source scenarios. Attached Figure Description

[0053] Figure 1 This is a schematic diagram of a host sequence filtering method for phage metagenomic sequencing data according to this application;

[0054] Figure 2 This is a schematic diagram of a host sequence filtering system for phage metagenomic sequencing data according to this application. Detailed Implementation

[0055] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0056] Example 1

[0057] Please see Figure 1 As shown in this embodiment, a method for filtering host sequences in phage metagenomic sequencing data includes:

[0058] S1. Obtain the phage metagenomic sequencing data of the sample to be analyzed and preprocess it to obtain quality filtered read data; extract the compositional features of the quality filtered read data and calculate the corresponding host abundance background value of the sample to be analyzed, while configuring the corresponding host genome library and setting the comparison judgment threshold.

[0059] S2. Map and align the quality-filtered read data to the host genome database, extract the read mapping status features of each read in the mapping and alignment results, perform source classification on the quality-filtered read data in combination with the alignment judgment threshold, and output the host mapped read set, the phage candidate read set, and the boundary candidate read set;

[0060] S3. Extract local breakpoints and pairing conflict features from the boundary candidate read set, segment the fragment sequence based on the local breakpoints, perform joint source re-determination on the segmentation results in combination with the pairing conflict features, identify the rescue fragment that belongs to the non-host sequence, and transfer the rescue fragment into the phage candidate read set;

[0061] S4. Assemble the sequence of the candidate reads of the transferred phage to obtain candidate contigs; determine the true source of the candidate contigs, establish a mapping link between the candidate contigs and the corresponding supporting reads, and make reverse corrections to the source classification of the corresponding supporting reads based on the determination results.

[0062] S5. Based on the reverse correction results, separate host-origin reads from the quality-filtered read data and output the filtered phage metagenomic sequencing data; calculate the read migration rate during the reverse correction process, and update and compensate the comparison threshold for subsequent batches of samples to be analyzed based on the deviation between the read migration rate and the corresponding host abundance background value.

[0063] In this embodiment, physical or biological samples are first collected from a pre-selected target environment. The samples are then enriched using physical methods such as centrifugation and ultrafiltration. Sequencing data of the corresponding samples is obtained using relevant equipment in the field of genetic engineering, resulting in phage metagenomic sequencing data. Sequence trimming algorithms are used to identify and remove sequencing adapters remaining at both ends of each original read in the phage metagenomic sequencing data. Low-quality bases are then screened out, and invalid reads that are too short after truncation or contain too many unidentified bases are removed. This preprocessing process yields higher-quality filtered read data.

[0064] Methods for extracting compositional features from quality-filtered reads and calculating the corresponding host abundance background value for the sample to be analyzed include:

[0065] Sequence sampling is performed on the quality filtered read data to identify the percentage of specific bases and fixed-length nucleotide short strings in each sampled read, which are then combined to form compositional features.

[0066] In this embodiment, sequence sampling refers to randomly selecting a portion of reads from the quality-filtered read data as sampled reads according to a preset sampling ratio; wherein the preset sampling ratio is set based on the magnitude of the sequencing data and the total computing resources to ensure that the quantity conforms to historical sampling experience.

[0067] In this embodiment, the percentage of a specific base content refers to the percentage of the sum of the number of guanine and cytosine bases in each sampled read relative to the total number of bases in the corresponding read. Different species’ genomes have relatively fixed distribution characteristics of guanine and cytosine content, while host bacteria and bacteriophages often have differences in the content of these bases. Therefore, the percentage of these two bases can be used as an indicator to distinguish different species origins.

[0068] Fixed-length nucleotide short strings refer to continuous subsequence fragments obtained by cutting the base sequence of each sampled read into a sliding window according to a preset fixed base length. The preset fixed base length is set based on practical experience in the field of bioinformatics to ensure that sufficient species discrimination is preserved while avoiding excessive length that would lead to a low matching hit rate.

[0069] A fixed-length nucleotide short string is input into a preset host feature hash index for matching, and the number of read segments that match is counted.

[0070] The preset host feature hash index library refers to a query database constructed based on the reference genome sequence of known host species. It obtains several short nucleotide strings by dividing the human genome, E. coli genome, or other known host bacterial genomes according to a certain fixed base length, and stores and indexes them in a hash table data structure.

[0071] By matching fixed-length nucleotide short strings with a preset host feature hash index table, when any one or more nucleotide short strings generated from the cutting of a sampled read segment match the preset host feature hash index, the sampled read segment is recorded as a match hit; the match hit counts of all sampled read segments are integrated to obtain the total number of read segments that have a match hit.

[0072] The initial ratio is obtained by calculating the ratio of the number of matched read segments to the total number of corresponding sampled read segments.

[0073] The initial proportion is obtained by calculating the ratio of the number of sampled reads that have a matching hit to the total number of all sampled reads. This value is used to indicate what proportion of reads in the current subset of sampled reads carry sequence features similar to those of the known host genome.

[0074] Based on the percentage of specific base content in each sampling segment, an overall base content distribution curve is plotted, the overlap area between the overall base content distribution curve and the corresponding base distribution interval is identified, and the base overlap is output.

[0075] In this embodiment, by using the percentage of a specific base content as the horizontal axis and the number of sampled reads with that percentage as the vertical axis, a total base content distribution curve is plotted for the percentage of a specific base content in all sampled reads. This curve is used to reflect the overall distribution of all reads in the current sample at the base composition level.

[0076] The corresponding base distribution interval refers to the percentage distribution range of the corresponding base content under standard conditions, determined based on the theoretical genome information of the known host. For example, the guanine and cytosine content in Escherichia coli is generally concentrated in the range of 50% to 52%, so this range is recorded as the corresponding base distribution interval.

[0077] By overlaying the theoretical base distribution range of the host with the actual overall base content distribution curve, and calculating the area enclosed by the overlapping part, the overlap area is obtained. This area is used as the base overlap degree. The larger the value, the more similar the base composition of the current sample is to the base composition of the known host.

[0078] The initial ratio was weighted and corrected using this base overlap to obtain the host abundance background value of the corresponding sample to be analyzed.

[0079] Because the initial ratio may contain conserved fragments leading to false positives or non-host sequences hitting short host sequences during the calculation process, resulting in an overestimation, the base overlap is used as a weighting factor and multiplied with the initial ratio to obtain the final host abundance background value. When the base overlap is high, it indicates that the macroscopic base composition conforms to the host characteristics, and the final host abundance background value is close to the initial ratio, indicating that the initial ratio estimation is reliable. When the base overlap is low, it indicates that the macroscopic base composition deviates from the standard base composition of the host, and the obtained host abundance background value is correspondingly reduced.

[0080] In this embodiment, the host abundance background value is used to represent the corrected estimate of the proportion of host-derived sequences in the overall sequencing data of the current sample to be analyzed, and characterizes the intensity of host contamination background that the corresponding sample to be analyzed has before formally entering the host filter.

[0081] It should be noted that the sampled reads obtained after performing sequence sampling on the quality-filtered read data in this embodiment are not used to describe the local attributes of a single read, but rather serve as a representative subset of the overall read set of the sample to be analyzed. Since the metagenomic sequencing data comes from the random fragmentation and random reading of sample nucleic acid fragments, under the same uniform sampling conditions, the sampled reads can statistically approximate the relative proportion of host sequences and non-host sequences and the base composition distribution in the overall read set. Therefore, the final host abundance background value can reflect the host background intensity of the overall sequencing data.

[0082] The methods for extracting the read segment mapping status features from the mapping comparison results and combining them with the comparison judgment threshold to perform source classification on the quality-filtered read segment data include:

[0083] The genome sequence in the host genome library is identified as a reference index. The quality-filtered read data is matched with the reference index to generate sequence alignment records as mapping comparison results.

[0084] In this embodiment, the host genome library refers to a database constructed based on existing theoretical knowledge that contains theoretical genome sequences of the corresponding species and records the theoretical genome attribute information of the host of the corresponding sample and other parameters related to the corresponding host genes. The complete genome sequence in the host genome library is used as a reference index. Each read in the quality-filtered read data is matched with the reference index to determine whether each read can be mapped to a certain position on the host genome sequence. The mapping and alignment results are recorded as sequence alignment records, which record the mapping position, orientation, base matching or mismatch status, and whether there are any deletions for each read on the host's theoretical genome.

[0085] The sequence coverage ratio, matching similarity, sequence alignment features with soft shearing markers, and position features of the paired read segments in the parsing mapping and alignment results are combined to form the read segment mapping state features.

[0086] The sequence coverage ratio refers to the proportion of the number of bases in the corresponding read that are successfully mapped to the host's theoretical genome out of the total number of bases in the read. This value is used to measure how much of the length of the corresponding read is covered by the host's theoretical genome sequence. The larger the value, the more completely the corresponding read falls on a certain region of the host genome.

[0087] Matching similarity refers to the proportion of bases in the successfully mapped base segment of the corresponding read segment that are completely identical to the corresponding bases in the host's theoretical genome, out of the total number of bases in the mapped segment. It is used to indicate the degree of similarity between the mapped part and the host's theoretical genome sequence at the base precision level. The larger the value, the stronger the local homology between the corresponding read segment and the host's theoretical genome sequence.

[0088] Sequence alignment features containing soft cut markers refer to the fact that some reads in the sequence alignment record do not fully participate in the mapping. Instead, there are base fragments at one or both ends of the read that are not included in the mapping comparison. These fragments are marked with soft cut markers in the sequence alignment record. The presence of this marker means that part of the base sequence of the corresponding read can be mapped to the host's theoretical genome, while the other part of the base sequence cannot be found at the corresponding position in the host's theoretical genome. This indicates that the read may cross the junction between the host sequence and the non-host sequence.

[0089] Paired-end read position characteristics refer to the sequence alignment record of paired reads generated by paired-end sequencing, which simultaneously contains the mapping position coordinates, orientation, and distance information of each end of the read. This distance information is the length of the inserted fragment between the two mapping positions. Under normal circumstances, the two ends of a paired read should be mapped to two positions on the host genome that are within one insertion fragment length apart in opposite directions. When there is a case where one end is mapped but the other end is not mapped, or when the mapping direction of both ends is abnormal or the length of the inserted fragment deviates significantly from the normal range, it is considered an abnormality of paired-end read position characteristics.

[0090] Read segments with both sequence coverage ratio and matching similarity not lower than the corresponding alignment judgment threshold are added to the host mapping read segment set.

[0091] The alignment threshold refers to an adaptive parameter configured based on the host abundance background value, which is consistent with the corresponding host. It includes two thresholds set for sequence coverage ratio and matching similarity. If the host abundance background value is large, it means that the background intensity of the host of the current sample is large. Therefore, the corresponding alignment threshold is relaxed to accelerate the host filtering process. Conversely, the alignment threshold is appropriately tightened to reduce damage to the phage sequence. The specific adjustment values ​​and settings are based on existing theoretical knowledge of the relevant sample species and hosts, and conform to the experience of historical host filtering. If the sequence coverage ratio and matching similarity of a certain read are not lower than the corresponding threshold, it means that the read is highly consistent with the theoretical genome of the host in terms of sequence coverage length and base matching accuracy. Therefore, it can be determined to be a sequence originating from the host.

[0092] Reads that could not be mapped to the host genome library at either end were added to the phage candidate read set.

[0093] For paired reads obtained from paired-end sequencing, if neither end of the read pair has any valid mapping position on the host's theoretical genome, it indicates that the base sequence of the read is not homologous to the known host's theoretical genome. It should be noted that if the sequence coverage ratio and matching similarity are lower than the corresponding alignment judgment threshold, the corresponding read should also be judged as unmapped. For unmapped reads, it is initially estimated that they are more likely to come from non-host species, so such reads are first added to the phage candidate read set.

[0094] Unmapped reads with soft shearing markers in the sequence alignment features, and reads with one end mapped and the other end unmapped or of abnormal length in the position features of the paired reads, are added to the boundary candidate read set.

[0095] In this embodiment, the boundary candidate read set refers to the set of reads that cannot be directly determined to be of host origin or non-host origin, representing a fuzzy classification boundary. If a soft splicing marker is present, it indicates that a portion of the base sequence of such paired-end reads may have successfully mapped to the host's theoretical genome, while the remaining base sequence did not participate in the mapping and may carry phage sequence information. If only one end of the paired-end reads is successfully mapped, it is considered a single-end mapping state. If the length of the inserted fragment between the two mapped positions significantly exceeds the standard inserted fragment length range under theoretical conditions, it is considered a length abnormality state. Reads exhibiting either of these abnormalities are added to the boundary candidate read set.

[0096] Methods for segmenting sequences based on local breakpoints include:

[0097] Traverse the sequence alignment records corresponding to each read segment in the boundary candidate read set and identify the soft cut markers within them.

[0098] The process involves reading the sequence marker of each read segment in the boundary candidate read segment set one by one, and then using the sequence marker as an index to query whether the corresponding sequence entry in the sequence alignment record has a soft cut marker.

[0099] Simultaneously, the alignment direction information and inserted segment length information of the double-ended read segments are extracted and combined as pairing conflict features.

[0100] The alignment direction information refers to the alignment direction of each end of the paired read when it is mapped to the host theoretical genome. Under normal circumstances, the two ends of the paired read should be mapped to the host theoretical genome in opposite directions, with one end being a forward mapping and the other end being a reverse mapping. The insertion fragment length information refers to the distance between the two ends of the paired read and their mapped positions on the host theoretical genome. The alignment direction information and the insertion fragment length information are combined into a pairing conflict feature, which serves as an abnormal feature indicator of paired reads.

[0101] Identify the cleaved base sequence indicated by the soft shearing marker and determine whether the number of bases in the cleaved base sequence has reached the preset lower limit of the truncation length.

[0102] In this embodiment, the cleaved base sequence refers to the original base sequence fragment extracted from the sequence alignment record that did not participate in host alignment within the base segment marked by the soft cleavage marker; the preset truncation length lower limit refers to the minimum effective number of bases threshold used to filter extremely short soft cleavage fragments caused by sequencing errors or random noise. The specific value is set based on filtering experience in the field of bioinformatics, for example, it is set to 30 bases or 50 bases.

[0103] If the number of bases reaches the preset truncation length limit, the specific boundary coordinates of the corresponding read segment will be identified as a local breakpoint based on the sequence alignment record.

[0104] If the number of bases reaches the preset truncation length lower limit, it indicates that the corresponding cleaved base sequence has sufficient sequence information to support subsequent operations. Conversely, if the number of bases is insufficient, it indicates that the corresponding sequence has insufficient information, which is likely caused by errors or noise and is not suitable for subsequent operations. The specific junction coordinates refer to the boundary positions between the base segments in the cleaved base sequence that are successfully mapped to the host theoretical genome and the unmapped base segments marked by soft cut markers. By scanning the corresponding read segments along a fixed direction according to the order of each read segment in the sequence alignment record, when the corresponding read segment has the first base position that changes from the mapped state to the soft cut state, this position is taken as the specific junction coordinates and confirmed as a local breakpoint.

[0105] Based on local breakpoints, the main body of the corresponding read segment is physically truncated, the unmapped short sequences are separated, and they are output as independent segments.

[0106] In this embodiment, physical truncation refers to separating the complete base sequence of the corresponding read segment at the location of the local breakpoint as the dividing line. One side is the successfully mapped host matching segment, and the other side belongs to the unmapped segment marked by the soft shearing mark, making it a short base sequence that can independently participate in subsequent operations, i.e., an independent segmentation fragment.

[0107] Methods for performing joint source re-determination on the segmentation results include:

[0108] Identify the paired end read segments corresponding to the independent segmentation segments, and obtain the paired end alignment direction information and paired end inserted segment length information from the paired end conflict features corresponding to the paired end read segments.

[0109] Since paired reads in paired-end sequencing share the same or associated sequence marker prefixes, it is possible to search for the other end read that forms a pair with the original source region in the quality-filtered read data by parsing the sequence markers inherited by the independently segmented fragments. At the same time, the alignment direction information and the length information of the inserted fragment at the paired end of the paired end read are obtained simultaneously as the data basis for subsequent operations.

[0110] When the paired end read is not mapped to the host genome library, and the paired end alignment direction information is reversed or the length information of the paired end inserted fragment exceeds the preset normal length range, the corresponding independent segment and the paired end read will be combined into a joint target.

[0111] The process begins by detecting whether the paired end reads are in an unmapped state. Simultaneously, it identifies the alignment direction and inserted fragment length information of the paired end. If the direction is reversed, it indicates that the relative mapping direction of the paired reads does not conform to the normal opposite-direction mapping relationship, suggesting that an external sequence may have been inserted in this region, causing the direction to flip. If the length of the inserted fragment exceeds the preset normal length range based on theoretical conditions, it indicates that the two ends may be separated by an external sequence not belonging to the host. If, while the paired end reads are in an unmapped state, the alignment direction is reversed or the inserted fragment length exceeds the limit, the corresponding independent segment and the paired end reads can be combined into an indivisible analysis unit, i.e., a joint target, serving as one of the data bases for subsequent operations.

[0112] The sequence composition of the joint target is extracted, and the corresponding region of the sequence composition is compared with the pre-defined conserved region of phage features for homology.

[0113] The process involves extracting the base sequences of the independently segmented fragments and the paired end reads contained in the joint target, and splicing or arranging them according to their relative positions on the original genome to form a base region that can represent the sequence information of the corresponding joint target. The pre-defined phage characteristic conserved region refers to a set of characteristic nucleotide sequence fragments with high sequence conservation in the known phage genome. That is, a sequence segment with high sequence similarity, low variation frequency and stable function in multiple phage genomes. Compared with other segments, the conserved region is more likely to maintain relatively consistent nucleotide arrangement characteristics among different phage species or strains.

[0114] Homology alignment refers to the local sequence alignment of the sequence composition of the joint target with the pre-defined conserved regions of phage characteristics to detect whether there are matching fragments in the base sequence of the joint target that have obvious homology with the pre-defined conserved regions of phage characteristics.

[0115] When a joint target contains a pre-defined conserved region of phage features, the joint target is identified as a salvage fragment of a non-host sequence, and the salvage fragment is transferred as a whole to the phage candidate read set.

[0116] If the homology alignment results show that the sequence composition of the joint target contains a matching fragment that is homologous to the pre-defined conserved phage characteristic region, it indicates that the joint target carries the characteristic gene information of the phage at the sequence level, and it can be confirmed that it originates from the phage rather than the host. At this time, the joint target as a whole is judged as a rescue fragment of the non-host sequence. The rescue fragment refers to the sequence that was classified into the boundary candidate read set in the original initial classification process because it contained part of the host matching segment. Through joint source re-determination, the non-host part of it is successfully identified and separated from the host matching segment and added to the phage candidate read set.

[0117] Methods for determining the true source of candidate contigs include:

[0118] The transferred salvage fragments are spliced ​​together with the original set of phage candidate reads into a predefined metagenomic assembly structure to generate a candidate contiguous group with a continuous sequence structure.

[0119] The preset metagenomic assembly structure refers to the computational module used to perform sequence assembly. This module extends and splices a large number of short reads from the input salvage fragment and the original phage candidate read set according to the overlap relationship between the base sequences, and merges short reads with first and last overlap into longer continuous sequences step by step to obtain candidate contigs.

[0120] Gene coding sequences extracted from candidate contigs are compared with host genome libraries for long-scale homology.

[0121] In this embodiment, the gene coding sequence refers to the potential protein-coding region with complete start and stop codons identified by the open reading frame prediction algorithm on the continuous base sequence of the candidate contiguous group. Since the full-length sequence of the candidate contiguous group may include both coding and non-coding regions, the coding region carries more specific species biological information and is more suitable as a basis for determining the origin. The corresponding overall candidate contiguous group is used as the query sequence and sequence alignment is performed with the host genome library at a larger sequence scale. Compared with the alignment at the short read level, the host homology signals that were originally scattered in multiple short reads are spliced ​​into continuous long fragments, making it easier to reveal their overall homology with the host genome.

[0122] The proportion of host homologous sequences corresponding to the gene coding sequence is statistically analyzed, and the sequence complexity and viral motif density of unmapped regions in the corresponding candidate contiguous groups are obtained.

[0123] The proportion of host homologous sequences refers to the ratio of the length of the gene coding sequence extracted from the candidate contiguous group to the total length of the gene coding sequence in which a clear homologous base match can be found in the host genome library; the larger this value is, the more components from the host are in the coding region of the candidate contiguous group.

[0124] Sequence complexity refers to the degree of diversity and randomness of base arrangements in unmapped regions. By statistically analyzing the frequency distribution of different base combination patterns in unmapped regions, a high degree of repetition or a single pattern in base arrangements indicates low sequence complexity; while a high degree of diversity in base arrangements without obvious repetition patterns indicates high sequence complexity.

[0125] Viral signature motif density refers to the frequency of occurrence of motifs in unmapped regions that match known short conserved nucleotide patterns specific to the phage genome. These motifs include, for example, characteristic short sequences that appear repeatedly in the phage genome but are relatively rare in the host genome, such as phage genome terminal repeat patterns, phage packaging signal sequence patterns, and phage promoter recognition sequence patterns. The higher the density, the higher the concentration of phage signatures carried in the unmapped regions, and the higher the probability that the corresponding candidate contiguous groups originate from real phages.

[0126] When the proportion of host homologous sequences is higher than the preset contamination threshold, and the sequence complexity of the unmapped region is lower than the expected complexity threshold, the corresponding candidate contiguous group will be identified as a host contamination contiguous group.

[0127] The preset contamination threshold is a critical value used to determine whether a candidate contiguous group is of host origin. The specific value is set based on historical comparison experience between known host genomes and known bacteriophage genomes. If the proportion of host homologous sequences is higher than the preset contamination threshold, it indicates that a large proportion of segments in the coding region of the corresponding candidate contiguous group have high homology with the host genome, and the probability of originating from the host is high.

[0128] The expected complexity threshold is a critical value used to help determine whether an unmapped region belongs to a low-complexity repetitive sequence in the host genome. It is set based on the complexity distribution of repetitive and low-complexity segments in the known host genome. Low sequence complexity indicates that the corresponding unmapped region is more likely to belong to a low-complexity segment in the host genome that failed to be mapped during the short read alignment stage due to excessive repetition.

[0129] When both threshold comparison conditions are met simultaneously, it indicates that the coding genes are largely homologous to the host, and the non-coding genes are more like repetitive sequences of the host than bacteriophages. Therefore, the corresponding candidate contigs are identified as host contamination contigs.

[0130] When the proportion of host homologous sequences is lower than the safety threshold and the density of viral characteristic motifs is higher than the significant density, the corresponding candidate contigs are determined to be true phage contigs.

[0131] The safety threshold is a critical value used to determine whether a candidate contiguous group can be safely identified as a host homologous sequence from a non-host source. This value is much smaller than the preset contamination threshold, and the specific value is set based on the distribution of a small number of conserved genes shared with the host in the phage genome in historical records.

[0132] Significant density refers to the density threshold used to determine whether the phage features in the unmapped region have reached a credible level. It is set based on the typical distribution of characteristic motifs in the known phage genome. If the proportion of host homologous sequences is lower than the safety threshold, it indicates that the corresponding candidate contiguous group has weak homology with the host genome and cannot support the judgment of host origin.

[0133] If the viral feature motif density is higher than the significant density, it indicates that the phage's unique sequence pattern is densely distributed in the unmapped region. Therefore, when both threshold comparison conditions are met at the same time, it means that the corresponding candidate contig is neither like the host structure nor contains the characteristics of the phage. Therefore, the candidate contig is determined to be a real phage contig.

[0134] The methods for establishing a mapping link between candidate contiguous groups and corresponding supporting reads, and for making reverse corrections to the source classification of the corresponding supporting reads based on the judgment results, include:

[0135] Extract the set of sequence identifiers for all initial short reads that construct each candidate contiguous group, and use each initial short read as a supporting read.

[0136] The initial short reads refer to the original short reads that actually participate in the splicing process and ultimately form the candidate contiguous groups during the previous assembly process; the sequence identifier set refers to the set of unique sequence markers carried by each initial short read; and by using each initial short read as a support read, it is shown that these short reads are the underlying building blocks that constitute the corresponding candidate contiguous groups.

[0137] Based on the affiliation relationship between the sequence identifier set and the corresponding candidate overlapping group, a mapping link is constructed from the candidate overlapping group to each corresponding supporting read segment in the corresponding sequence identifier set.

[0138] The attribution relationship refers to the subordinate correspondence between the sequence identifier of each supporting read segment and the candidate contiguous group it constructs. That is, one candidate contiguous group corresponds to multiple supporting read segments, and one supporting read segment belongs to one or more candidate contiguous groups. Based on the attribution relationship, a one-way retrieval link is established from the candidate contiguous group to its underlying supporting read segment, that is, the mapping link.

[0139] When a candidate contiguous group is identified as a host-contaminated contiguous group, a reverse tracing operation is performed on the corresponding mapping link to lock all abnormal support reads that constitute the host-contaminated contiguous group.

[0140] The reverse source tracing operation refers to extracting all supporting reads from the sequence identifier set corresponding to the candidate contigs identified as host contamination contigs based on the candidate contigs, through the corresponding mapping links. Abnormal supporting reads are those that were classified into the phage candidate read set during the historical screening process, but after long-sequence large-scale homology comparison, the candidate contigs formed were confirmed as the source of host contamination. This indicates that the source classification of these supporting reads is incorrect and they should actually belong to the host sequence. Because the information content of short reads is limited, they were not accurately screened during the screening process.

[0141] The source classification labels of all abnormal support reads were uniformly changed to host classification, and all abnormal support reads were separated from the phage candidate read set after being transferred to the salvage fragment and merged and transferred to the host-mapped read set.

[0142] The source classification label refers to the source attribute marker assigned to each read during the initial screening process, i.e., non-host classification, and is assigned to the phage candidate read set or the boundary candidate read set. Since the abnormal supporting reads are found to actually belong to the host sequence after comparison, the source classification label is changed to the host classification. At the same time as modifying the label, these abnormal supporting reads are physically removed from the current phage candidate read set and added to the host-mapped read set.

[0143] The methods for calculating read migration rate during the reverse correction process, based on the deviation between read migration rate and corresponding host abundance background value, include:

[0144] The total number of anomalous support reads that are separated from the phage candidate read set after being transferred to the rescue fragment and transferred to the host-mapped read set during the statistical reverse correction process.

[0145] The total number of sequences refers to the cumulative number of all reads that were forcibly separated from the phage candidate read set and finally merged into the host mapped read set during the reverse correction process through reverse tracing of the mapping link; this value directly reflects how many reads that originally belonged to the host were not correctly classified in the initial screening stage.

[0146] The read migration rate is obtained by calculating the ratio of the total number of reads in the quality filtered read data to the total number of reads in the sequence. The read migration rate is used to convert the scale of read migration during the reverse correction process into a proportion relative to the total amount of data, which serves as the data basis for subsequent operations.

[0147] Extract the host abundance background value of the corresponding sample to be analyzed, and obtain the preset escape tolerance coefficient. Calculate the product of the preset escape tolerance coefficient and the host abundance background value to obtain the maximum allowable tolerance value of the corresponding sample to be analyzed.

[0148] In this embodiment, the preset escape tolerance coefficient is a proportional coefficient used to convert the host abundance background value into an upper limit of the allowable missed filtering ratio. Since the information content of short reads is limited and the host genome library may not be completely covered in the actual host sequence filtering process, a small number of host reads may pass the initial screening and be corrected back in the reverse correction. Therefore, the value of this proportional coefficient cannot be too large, otherwise even if a large number of host sequences escape, no alarm will be triggered, and the subsequent threshold will never be optimized. At the same time, the threshold cannot be set too small, otherwise even a small amount of missed filtering within the normal range will frequently trigger threshold adjustment, resulting in unstable fluctuations in parameters. Therefore, the specific value of this coefficient is set based on the distribution of the missed filtering probability of the same type of sample in the historical record.

[0149] The maximum allowable tolerance value refers to the upper limit of read migration rate that can be tolerated under the current host background intensity of the sample to be analyzed. It is obtained by calculating the product of the host abundance background value and the preset escape tolerance coefficient. For samples with higher host abundance background values, the cardinality of the host sequence in the overall data is larger. Even if the alignment judgment threshold is set reasonably, it is easier to generate more missed reads due to the excessive cardinality. Therefore, a larger tolerance value space is required. Conversely, for samples with lower host abundance background values, the cardinality of the host sequence is smaller. If a high read migration rate still occurs, it further indicates that the alignment judgment threshold applied in the initial screening is inappropriate. Therefore, the tolerance value space should be smaller. The tolerance value is synchronously scaled with the host background value through product calculation.

[0150] If the read segment migration rate is higher than the maximum allowable tolerance value, the corresponding deviation state is determined to be a positive deviation state; otherwise, the deviation state is determined to be a stable deviation state.

[0151] The deviation state refers to the quantitative assessment conclusion of the rationality of the alignment judgment threshold corresponding to the preliminary screening of the current sample to be analyzed. If the read migration rate is higher than the maximum allowable tolerance value, it means that the proportion of reads corrected back to the host source in the reverse correction has exceeded the normal tolerance range of the current host background intensity. This indicates that the alignment judgment threshold used in the preliminary screening is set too loosely, resulting in a large number of host sequences not being effectively intercepted in this screening stage. Therefore, this deviation state is judged as a positive deviation state, indicating that the alignment judgment threshold of the batch corresponding to the current sample to be analyzed needs to be reduced in subsequent batches.

[0152] When the read migration rate is not higher than the maximum allowable tolerance value, it means that the proportion of reads corrected in the reverse correction is within the normal tolerance range of the current host background intensity. This indicates that the comparison judgment threshold used in the initial screening has a basically reasonable interception strength of the host sequence, and there is no need to make significant adjustments to the threshold in the future. Therefore, this deviation state is judged as a stable deviation state.

[0153] The methods for updating and compensating for the comparison threshold of subsequent batches of samples to be analyzed include:

[0154] If the deviation is positive, obtain the sample source label of the corresponding sample to be analyzed, and calculate the ratio of the difference between the corresponding read segment mobility and the maximum allowable tolerance value.

[0155] The sample source label refers to the label used to mark the source of the sample to be analyzed, such as human intestinal sample, animal intestinal sample, and sewage environmental sample. In this embodiment, the difference ratio refers to the proportion of the extra part of the corresponding read migration rate that exceeds the maximum allowable tolerance value to the maximum allowable tolerance value. The larger the value, the further the actual leakage level deviates from the normal range, and therefore the greater the compensation required.

[0156] Calculate the product of the difference ratio and the preset adjustment step size to obtain the compensation value.

[0157] The preset adjustment step size refers to the proportional parameter used to control the magnitude of a single threshold compensation. The specific value is set based on the threshold adjustment range of similar samples in the historical record. The adjustment magnitude cannot be too large or too small. A large adjustment magnitude can easily cause the threshold of subsequent batches to be adjusted from one extreme to another, resulting in instability in the screening process. A small adjustment magnitude may require the accumulation of multiple batches to correct the threshold to a normal level. By calculating the product of the difference ratio and the preset adjustment step size, the specific correction amount that needs to be superimposed on the original comparison judgment threshold in this adjustment is obtained, so that the final compensation value is constrained by both the severity of the deviation and the adjustment step size.

[0158] The corresponding comparison judgment threshold and the compensation value are superimposed to obtain the compensated threshold. The calculated compensation value is then arithmetically added to the original comparison judgment threshold to obtain the updated compensated threshold.

[0159] The compensated threshold is then bound to the corresponding sample source label and saved.

[0160] By binding the compensated threshold to the corresponding sample source label, it is ensured that the compensation result of the threshold can only be effective for subsequent samples with the same biological source type. For example, if a positive deviation is found during the processing of a batch of human intestinal samples and a compensated threshold is calculated, this threshold can only be used for subsequent human intestinal samples and cannot be applied to other types of samples such as sewage or soil samples, because samples from different sources have completely different host composition structures and host abundance characteristics.

[0161] When subsequent batches with the same sample source label appear, the corresponding compensated threshold is used as the comparison judgment threshold.

[0162] When a new sample to be analyzed is received, the sample source label of the sample is read first. If the label can be matched, the corresponding compensated threshold is extracted to replace the initially generated comparison judgment threshold. If it cannot be matched, the comparison judgment threshold of the previous batch is used. This enables the entire host sequence filtering process to have cross-batch adaptive learning capability for samples of the same type, improving the efficiency and accuracy of the overall process.

[0163] This embodiment implements a host sequence filtering method for phage metagenomic sequencing data by acquiring phage metagenomic sequencing data of the samples to be analyzed and performing host sequence filtering. Compared with the prior art, by specifically configuring the host genome library and the alignment judgment threshold, it improves the situation where the existing technology uses fixed filtering conditions for samples with different host backgrounds, which leads to mismatch problems. By performing fragment sequence segmentation, it avoids the overall deletion of boundary reads that contain both host and non-host regions, preventing the loss of effective phage fragments. By performing joint source re-determination, it solves the problem of broken evidence chains of structural anomalies caused by splitting paired reads for independent processing. By performing true source determination and reverse correction, it enables host contamination sequences missed in the initial screening stage to be re-identified and recovered at the long sequence scale. By updating and compensating the alignment judgment threshold of subsequent batches based on the deviation between read migration rate and host abundance background value, it improves the adaptability of the host sequence filtering process in cross-batch applications of the same type. It improves the accuracy of host sequence filtering for phage metagenomic sequencing data, and enhances the stability and generalization ability of the filtering results under different sample source scenarios.

[0164] Example 2

[0165] Please see Figure 2 As shown, for parts not described in detail in this embodiment, please refer to the description in Embodiment 1. A host sequence filtering system for phage metagenomic sequencing data is provided, comprising:

[0166] The preliminary filtering module acquires and preprocesses the phage metagenomic sequencing data of the sample to be analyzed to obtain quality-filtered read data; it extracts the compositional features of the quality-filtered read data and calculates the corresponding host abundance background value of the sample to be analyzed, while configuring the corresponding host genome library and setting the alignment judgment threshold.

[0167] The source classification module maps and aligns the quality-filtered read data to the host genome database, extracts the read mapping status features of each read in the mapping and alignment results, performs source classification on the quality-filtered read data in combination with the alignment judgment threshold, and outputs the host mapped read set, the phage candidate read set, and the boundary candidate read set.

[0168] The fragment transfer module extracts local breakpoints and pairing conflict features from the boundary candidate read set, performs fragment sequence segmentation based on local breakpoints, and performs joint source re-determination on the segmentation results in combination with pairing conflict features to identify rescue fragments belonging to non-host sequences, and transfers the rescue fragments into the phage candidate read set;

[0169] The reverse correction module assembles the sequence of the transferred phage candidate read set to obtain candidate contigs; it determines the true source of the candidate contigs, establishes a mapping link between the candidate contigs and the corresponding supporting reads, and performs reverse correction on the source classification of the corresponding supporting reads based on the determination results.

[0170] The compensation update module separates host-origin reads from the quality-filtered read data based on the reverse correction results and outputs the filtered phage metagenomic sequencing data; it calculates the read migration rate during the reverse correction process and updates and compensates the comparison threshold for subsequent batches of samples to be analyzed based on the deviation of the read migration rate from the corresponding host abundance background value; the modules are connected to each other via wired and / or wireless means.

[0171] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

[0172] All formulas in this manual are dimensionless and calculated numerically. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

[0173] Although embodiments of this application have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of this application, the scope of which is defined by the claims and their equivalents.

Claims

1. A method for filtering host sequences in phage metagenomic sequencing data, characterized in that, include: S1. Obtain the phage metagenomic sequencing data of the sample to be analyzed and perform preprocessing to obtain quality-filtered read data; Extract the compositional features of the quality filtered read data and calculate the host abundance background value of the corresponding sample to be analyzed. At the same time, configure the corresponding host genome library and set the comparison judgment threshold. S2. Map and align the quality-filtered read data to the host genome database, extract the read mapping status features of each read in the mapping and alignment results, perform source classification on the quality-filtered read data in combination with the alignment judgment threshold, and output the host mapped read set, the phage candidate read set, and the boundary candidate read set; S3. Extract local breakpoints and pairing conflict features from the boundary candidate read set, segment the fragment sequence based on the local breakpoints, perform joint source re-determination on the segmentation results in combination with the pairing conflict features, identify the rescue fragment that belongs to the non-host sequence, and transfer the rescue fragment into the phage candidate read set; S4. The candidate read set of the transferred phage is sequence assembled to obtain the candidate contiguous group; The candidate contiguous groups are determined to have the true source, a mapping link between the candidate contiguous groups and the corresponding supporting reads is established, and the source classification of the corresponding supporting reads is reversed based on the determination results. S5. Based on the reverse correction results, separate the host-origin reads from the quality-filtered read data and output the filtered phage metagenomic sequencing data; The read migration rate during the statistical reverse correction process is used to update and compensate the comparison threshold for subsequent batches of samples to be analyzed based on the deviation between the read migration rate and the corresponding host abundance background value.

2. The method for filtering host sequences in phage metagenomic sequencing data according to claim 1, characterized in that, The methods for extracting the compositional features of the quality filtered read data and calculating the corresponding host abundance background value of the sample to be analyzed include: Sequence sampling is performed on the quality filtered read data to identify the percentage of specific bases and fixed-length nucleotide short strings in each sampled read, which are then combined to form compositional features. A fixed-length nucleotide short string is input into a preset host feature hash index for matching, and the number of read segments that match is counted. The initial ratio is obtained by calculating the ratio of the number of matched reads to the total number of corresponding sampled reads; the overall base content distribution curve is plotted based on the specific base content percentage of each sampled read, the overlap area between the overall base content distribution curve and the corresponding base distribution interval is identified, and the base overlap is output; the initial ratio is weighted and corrected using the base overlap to obtain the host abundance background value of the corresponding sample to be analyzed.

3. The method for filtering host sequences in phage metagenomic sequencing data according to claim 2, characterized in that, The method of extracting the reading segment mapping status features of each reading segment in the mapping comparison results and performing source classification on the quality-filtered reading segment data in combination with the comparison judgment threshold includes: The genome sequence in the host genome library is identified as a reference index. The quality filtered read data is matched with the reference index to generate sequence alignment records as mapping comparison results. The sequence coverage ratio, matching similarity, sequence alignment features with soft shearing markers, and position features of the paired read segments in the parsing mapping and alignment results are combined into read segment mapping state features. Reads with sequence coverage and matching similarity not lower than the corresponding alignment threshold are added to the host-mapped read set; reads that cannot be mapped to the host genome library at either end are added to the phage candidate read set. Unmapped reads with soft shearing markers in the sequence alignment features, and reads with one end mapped and the other end unmapped or of abnormal length in the position features of the paired reads, are added to the boundary candidate read set.

4. The method for filtering host sequences in phage metagenomic sequencing data according to claim 3, characterized in that, The method of segmenting a sequence based on local breakpoints includes: Traverse the sequence alignment records corresponding to each read in the boundary candidate read set and identify the soft cut markers; at the same time, extract the alignment direction information and insertion fragment length information of the paired reads and combine them as pairing conflict features. Identify the cleaved base sequence indicated by the soft shearing marker and determine whether the number of bases in the cleaved base sequence reaches the preset truncation length lower limit; if the number of bases reaches the preset truncation length lower limit, then confirm the specific boundary coordinates of the corresponding read segment as a local breakpoint based on the sequence alignment record; Based on local breakpoints, the main body of the corresponding read segment is physically truncated, the unmapped short sequences are separated, and they are output as independent segments.

5. The method for filtering host sequences in phage metagenomic sequencing data according to claim 4, characterized in that, The method for performing joint source re-determination on the segmentation results includes: Identify the paired end read segments corresponding to the independent segmentation segments, and obtain the paired end alignment direction information and paired end inserted segment length information from the paired end conflict features corresponding to the paired end read segments; When the paired end read is not mapped to the host genome library, and the paired end alignment direction information is reversed or the length information of the paired end inserted fragment exceeds the preset normal length range, the corresponding independent segment and the paired end read will be combined into a joint target. The sequence composition of the joint target is extracted, and the corresponding region of the sequence composition is compared with the pre-defined conserved region of phage features. When the joint target contains the pre-defined conserved region of phage features, the joint target is identified as a rescue fragment of a non-host sequence, and the rescue fragment is transferred as a whole to the phage candidate read set.

6. The method for filtering host sequences in phage metagenomic sequencing data according to claim 5, characterized in that, The method for determining the true source of candidate overlapping groups includes: The transferred salvage fragments and the original set of phage candidate reads are input into a pre-defined metagenomic assembly structure and spliced ​​to generate a candidate contiguous group with a continuous sequence structure. Gene coding sequences extracted from candidate contigs are compared with host genome databases for long-scale homology. The proportion of host homologous sequences corresponding to the gene coding sequence is calculated, and the sequence complexity and viral motif density of unmapped regions in the corresponding candidate contigs are obtained. When the proportion of host homologous sequences is higher than the preset contamination judgment threshold, and the sequence complexity of the unmapped region is lower than the expected complexity threshold, the corresponding candidate contiguous group will be judged as the host contamination contiguous group. When the proportion of host homologous sequences is lower than the safety threshold and the density of viral characteristic motifs is higher than the significant density, the corresponding candidate contigs are determined to be true phage contigs.

7. The method for filtering host sequences in phage metagenomic sequencing data according to claim 6, characterized in that, The method of establishing a mapping link between candidate overlap groups and corresponding supporting reads, and then reversing the source classification of the corresponding supporting reads based on the determination results, includes: Extract the sequence identifier set of all initial short reads used to construct each candidate contiguous group, and use each initial short read as a supporting read; Based on the affiliation relationship between the sequence identifier set and the corresponding candidate overlapping group, a mapping link is constructed from the candidate overlapping group to each corresponding supporting read segment in the corresponding sequence identifier set; When the corresponding candidate contiguous group is determined to be a host-contaminated contiguous group, a reverse tracing operation is performed on the corresponding mapping link to lock all abnormal support reads that constitute the host-contaminated contiguous group. The source classification labels of all abnormal support reads were uniformly changed to host classification, and all abnormal support reads were separated from the phage candidate read set after being transferred to the salvage fragment and merged and transferred to the host-mapped read set.

8. The method for filtering host sequences in phage metagenomic sequencing data according to claim 7, characterized in that, The read migration rate in the statistical reverse correction process, based on the deviation between the read migration rate and the corresponding host abundance background value, includes the following methods: The total number of anomalous supporting reads that are separated from the phage candidate read set after being transferred to the rescue fragment and transferred to the host-mapped read set during the statistical reverse correction process; The read mobility is obtained by calculating the ratio of the total number of reads in the sequence to the total number of reads in the quality-filtered read data. Extract the host abundance background value of the corresponding sample to be analyzed, and obtain the preset escape tolerance coefficient. Calculate the product of the preset escape tolerance coefficient and the host abundance background value to obtain the maximum allowable tolerance value of the corresponding sample to be analyzed. If the read migration rate is higher than the maximum allowable tolerance value, the corresponding deviation state is determined to be a positive deviation state; otherwise, the deviation state is determined to be a stable deviation state.

9. The method for filtering host sequences in phage metagenomic sequencing data according to claim 8, characterized in that, The method for updating and compensating the comparison judgment threshold for subsequent batches of samples to be analyzed includes: If the deviation is positive, obtain the sample source label of the corresponding sample to be analyzed, and calculate the ratio of the difference between the corresponding read segment mobility and the maximum allowable tolerance value. Calculate the product of the difference ratio and the preset adjustment step size to obtain the compensation value; superimpose the corresponding comparison judgment threshold and the compensation value to obtain the compensated threshold; bind the compensated threshold with the corresponding sample source label and save it; when subsequent batches with the same sample source label appear, apply the corresponding compensated threshold as the comparison judgment threshold.

10. A host sequence filtering system for phage metagenomic sequencing data, used to implement the host sequence filtering method for phage metagenomic sequencing data according to any one of claims 1-9, characterized in that, include: The preliminary filtering module acquires and preprocesses the phage metagenomic sequencing data of the sample to be analyzed to obtain quality filtered read data; Extract the compositional features of the quality filtered read data and calculate the host abundance background value of the corresponding sample to be analyzed. At the same time, configure the corresponding host genome library and set the comparison judgment threshold. The source classification module maps and aligns the quality-filtered read data to the host genome database, extracts the read mapping status features of each read in the mapping and alignment results, performs source classification on the quality-filtered read data in combination with the alignment judgment threshold, and outputs the host mapped read set, the phage candidate read set, and the boundary candidate read set. The fragment transfer module extracts local breakpoints and pairing conflict features from the boundary candidate read set, performs fragment sequence segmentation based on local breakpoints, and performs joint source re-determination on the segmentation results in combination with pairing conflict features to identify rescue fragments belonging to non-host sequences, and transfers the rescue fragments into the phage candidate read set; The reverse correction module performs sequence assembly on the transferred set of phage candidate reads to obtain candidate contigs. The candidate contiguous groups are determined to have the true source, a mapping link between the candidate contiguous groups and the corresponding supporting reads is established, and the source classification of the corresponding supporting reads is reversed based on the determination results. The compensation update module separates host-origin reads from the quality-filtered read data based on the reverse correction results and outputs the filtered phage metagenomic sequencing data. The read migration rate during the statistical reverse correction process is used to update and compensate the comparison threshold of subsequent batches of samples to be analyzed based on the deviation of the read migration rate from the corresponding host abundance background value; the modules are connected to each other via wired and / or wireless means.

Images