Panel short query dual engine adaptive design method for capture sequencing probes for outgroup genomes

Abstract

The application provides a capture sequencing probe panel short query dual-engine adaptive design method for a close relative species genome. The method strictly defines the calculation rule of the judgment field, only calculates the percent identity and the query coverage based on the optimal local alignment segment itself, excludes the alignment block length as the coverage estimation basis from the source, avoids the systematic false positive caused by anchor extension in the short query scene based on the preset alignment tool of the short read, and realizes the bottom-up of the edge scene through the fallback engine, thereby significantly improving the reliability of the close relative species identification panel design.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics, and relates to the field of probe design for targeted capture sequencing, and more specifically to a short query dual-engine adaptive design method for capture sequencing probe panels for genomes of closely related species. Background Technology

[0002] Hybridization capture-based targeted sequencing (HMRS) uses biotin-tagged oligonucleotide probe panels (typically 80–120 nt in length) targeting specific regions of the genome to enrich DNA molecules containing the target sequence from the whole-genome background before sequencing. Compared to amplicon methods, HMRS offers advantages such as uniform coverage, good fragment length tolerance, and compatibility with low-quality fragmented DNA (e.g., cfDNA, FFPE samples), making it the mainstream approach for low-copy targets in molecular diagnostics (IVD) and clinical translation research.

[0003] Identifying closely related species is a typical application of capture sequencing. When a target species shares most of its genome sequence with multiple closely related species (such as different species within the same genus), a single species' whole-genome probe panel cannot determine its species affiliation. The panel must contain only species-specific regions of the target species relative to its closely related species; otherwise, the reaction will enrich sequences from closely related species simultaneously, making it impossible for downstream sequencing to distinguish between them.

[0004] For target species with a genome size of 1 Mb or more (such as parasite nuclear genomes of ~100 Mb, fungi of ~30 Mb, and bacteria of ~5 Mb), since commercial probe synthesis panels usually have an upper limit on total length, the total number of full-length tiling probes that completely cover the target genome will exceed the feasible scale for commercial synthesis. Furthermore, only some segments in a large genome are specific to closely related species. Therefore, panel design must first perform species-specific segment screening and then perform tiling on the selected segments.

[0005] The industry-standard species-specific region screening process typically includes: 1) cutting candidate window sequences from the target species' reference genome at fixed steps (e.g., every 100–500 nt); 2) aligning the candidate windows with the reference genomes of multiple closely related species to calculate identity and query coverage; 3) setting thresholds (typically, identity ≥ 90% and coverage ≥ 80% are considered shared, 70% ≤ identity < 90% are considered cross-reactivity, and the rest are considered species-specific); 4) retaining the windows determined to be species-specific as input target regions for subsequent probe tiling; 5) performing host (e.g., GRCh38) off-target checks on the target regions and removing repetitive element contamination to obtain the final panel.

[0006] However, alignment tools based on short-read presets are designed for high-throughput genome alignment of paired-end shortreads. Their internal parameters (including small k-mers, loose extension, and independence from pairing information) are designed for scenarios where the approximate mapping position is known and a 30-50 nucleotide consecutive perfect match is sufficient for anchoring. When the input is a single-end short window of ≤200 nt and the overall similarity between the query sequence and the reference sequence is in the 70–90% range (i.e., the typical identity range of closely related species), such presets will exhibit the following systematic bias (see details). Figure 2The algorithm selects locally perfectly matched segments of 30–50 nt as anchor points. When the anchor points are extended outward, the short read preset has a high tolerance for extended mismatches, and the extension cost function is biased towards extending the alignment block. The length of the final output alignment block is artificially lengthened (typically from 30 nt at the anchor point to more than 160 nt). The proportion of real matches in the alignment block is only 15% to 25%, but the alignment block-level identity field output by PAF (calculated by (block_len−NM) / block_len within the alignment block) appears to be 85% to 95%. The alignment block-level querycoverage field output by PAF is calculated as alignment block length / query length (instead of real matching segment length / query length). Therefore, a window that actually shares only about 30 nt will be reported as query coverage greater than or equal to 80%. The downstream threshold judgment (identity≥90% and coverage≥80% are considered shared) is deceived by this apparent high value, and the real species-specific window is misjudged as a shared window. Because many windows that should have been shared or cross-reactive were misjudged as shared and excluded, the set of untouched thresholds in the remaining window set (determined as species-specific according to complementary logic) actually contains a large number of false positive species-specific windows.

[0007] In addition, existing methods have the following problems: host reference genome off-target screening does not consider interference from nuclear mitochondrial transposon insertion fragments, leading to mitochondrial probes being misclassified as high-risk off-targets; target sources are singular, and multi-source information is not structurally fused; species-specific screening thresholds rely on single-point empirical values, and threshold sensitivity is not systematically scanned; high-copy repeat marker remapping is incomplete, only taking the best alignment result and losing copy information; probe-level host off-target classification granularity is insufficient; and probe-level GC content and homopolymer quality filtering are lacking. Summary of the Invention

[0008] To address the systematic false positives caused by anchor point extension in existing short-read-based alignment tools in short-query scenarios, as well as various systemic defects in panel design and related processes, this invention provides a dual-engine adaptive design method for short-query capture sequencing probe panels used in genomes of closely related species.

[0009] The technical solution adopted in this invention is: a short query design method for capture sequencing probe panels, including: S1. Obtaining the reference genome sequence of the target species and the reference genome sequence of closely related species; S2. Extracting candidate window sequences with a length not exceeding 200 nucleotides from the reference genome of the target species at a fixed step size; S3. Performing dual-engine alignment on the candidate window and the reference genome of closely related species, including:

[0010] S3a. The main engine based on the seed-expansion precise alignment algorithm is used for alignment. The first optimal local alignment segment that meets the preset e-value threshold is selected and the fields percent identity and query coverage are calculated, denoted as (pident_p, qcovs_p). The percent identity field is defined as the ratio of the number of matching bases in the optimal local alignment segment to the length of the optimal local alignment segment, and the query coverage field is defined as the ratio of the length of the optimal local alignment segment to the length of the candidate window sequence.

[0011] S3b. When any preset trigger condition is met, the fallback engine based on the assembly-level low difference preset is triggered to re-compare, select the second optimal local comparison segment and calculate the percent identity and query coverage according to the same field definition rules as the main engine, denoted as (pident_f, qcovs_f).

[0012] S3c. If sub-step S3a obtains the first optimal local alignment segment and does not meet any triggering condition of sub-step S3b, then (pident_p, qcovs_p) is used to execute step S4; otherwise, (pident_f, qcovs_f) of sub-step S3b is used to execute step S4.

[0013] S4. Record (pident_p, qcovs_p) or (pident_f, qcovs_f) output in step S3 as (identity, coverage), and make the following judgment for each candidate window based thereon: if the judgment for any closely related species reference genome satisfies identity≥T_shared and coverage≥T_cov, the fusion is determined as shared; otherwise, if the judgment for any closely related species reference genome satisfies T_cross≤identity<T_shared and coverage≥T_cov, the fusion is determined as cross_reactive; in other cases, the fusion is determined as species_specific; wherein T_shared, T_cross and T_cov are preset thresholds; S5. Collect candidate windows that are fusion-determined as species_specific in step S4 as targets, and perform tiling to generate a candidate oligonucleotide probe sequence set;

[0014] S6. Output the candidate oligonucleotide probe sequence set as a species-specific capture probe panel.

[0015] The invention provides a method for designing a capture sequencing probe panel for closely related species genomes. Under the scenario of short query window of ≤200 nt and cross-reference alignment of closely related species, by strictly defining the calculation rule of the judgment field, percent identity and query coverage are calculated only based on the optimal local alignment segment itself, which stably outputs an identity / coverage quantitative judgment that conforms to the real matching segment of the query window, eliminates the alignment block length as the estimation basis of coverage from the source, avoids systematic false positives caused by anchor extension in short query scenarios based on alignment tools preset for short reads, and realizes a fallback for edge scenarios through a fallback engine, which significantly improves the reliability of the design of the identification panel for closely related species. The output can be directly delivered to synthesize the species-specific capture probe sequence set, which is used for hybridization enrichment of downstream capture sequencing libraries.

[0016] Preferably, in step S1, the total sequence length of the target species reference genome is not less than 10 6 nt.

[0017] Preferably, in step S2, the fixed step size is 20~200 nt, and the candidate window sequence length is 50~200 nt, preferably 100~200 nt.

[0018] Preferably, in step S2, the candidate window sequence is subjected to low-complexity filtering before being cut out. The low-complexity filtering includes at least: dust filtering, tandem repeat filtering, and long homopolymer filtering.

[0019] Preferably, in step S3, the selection rule for the first optimal local alignment segment is: among all local alignment results output by the main engine on the candidate window and the reference genome pair of closely related species, the local alignment segment with the maximum bit score is taken as the first optimal local alignment segment.

[0020] As a preferred option, the selection rule for the second optimal local alignment segment is as follows: among all local alignment results output by the fallback engine on the candidate window and the reference genome pair of closely related species, the local alignment segment with the maximum value of the actual matching segment length is taken as the second optimal local alignment segment.

[0021] Preferably, the percent identity and query coverage fields are calculated based on the actual matching segments within the selected first optimal local alignment segment or second optimal local alignment segment, rather than using the total length of the alignment block as the basis for coverage estimation.

[0022] Preferably, in step S3, the main engine and the fallback engine are executed in parallel: the main engine and the fallback engine simultaneously compare the candidate window with the reference genome of the closely related species.

[0023] Preferably, in sub-step S3b, the preset triggering condition includes at least one of the following:

[0024] (i) The main engine did not output the first optimal local alignment segment;

[0025] (ii) The maximum bit score of the first optimal local alignment segment is lower than the preset lower limit of the bit score;

[0026] (iii) The maximum segment length of the first optimal local comparison segment is lower than the preset length lower limit;

[0027] (iv) The minimum e-value of the first optimal local alignment segment is greater than the preset e-value threshold;

[0028] The preset lower limit for bit score is 30-50, the preset lower limit for length is 20-30 nucleotides, and the preset e-value threshold is 10. -3 ~10 -8 .

[0029] Preferably, if neither the first optimal local alignment segment nor the second optimal local alignment segment is obtained in sub-steps S3a and S3b, then the candidate window is determined to have no significant homology with the reference genome of the closely related species.

[0030] Preferably, when the fallback engine still does not output a result that meets the preset threshold, the main engine is further invoked to re-align the candidate window with the reference genome of the closely related species using a shorter word length parameter (preferably -word_size 7 or -task blastn-short).

[0031] Preferably, in sub-step S3a, the main engine is selected from at least one of blastn, megablast, blastn-short, BLAT, LAST, USEARCH, edlib, and parasail.

[0032] Preferably, in sub-step S3b, the fallback engine is selected from at least one of minimap2 -cx asm5, -cxasm10, -cx asm20, LAST, and nucmer.

[0033] Preferably, in step S4, T_shared is 85%~95%, T_cross is 65%~80%, and T_cov is 70%~90%.

[0034] Preferably, in step S4, T_shared and T_cross are assigned values ​​based on the genome-wide average nucleotide similarity (ANI) between the target species and closely related species according to the following hierarchical mapping:

[0035] When ANI≥95%, T_shared takes a value of 95%~98%, and T_cross takes a value of 80%~90%.

[0036] When 80%≤ANI<95%, T_shared takes a value of 90%~95%, and T_cross takes a value of 70%~80%;

[0037] When 70%≤ANI<80%, T_shared takes a value of 85%~90%, and T_cross takes a value of 65%~75%;

[0038] When ANI < 70%, T_shared is set to 85%, T_cross is set to 65% (i.e., the upper limit of the default value), and T_cov is set to 65%~75% (preferably 70%) to avoid excessive underreporting;

[0039] When the ANI estimate falls within ±1% of the stratification boundary (e.g., 79%–81%, 94%–96%), the more stringent threshold combination between the two adjacent levels is used to avoid decision jumps caused by hard boundary shearing.

[0040] Preferably, after step S4 and before step S5, a threshold sensitivity scanning step is included: for the pre-target set composed of candidate windows determined to be species_specific, each pre-target is determined to have species specificity under a set of preset candidate thresholds, and a two-dimensional sensitivity matrix of target and threshold is output; a pre-target is included in the target of step S5 only when the determination results are consistent under all candidate thresholds. This invention ensures that the determination of species-specific targets in the final panel has cross-threshold robustness through threshold sensitivity scanning quality control.

[0041] Preferably, after step S4 and before step S5, a target region merging step is included: adjacent candidate window sequences that are both determined to be species_specific are merged into a single continuous target according to the following rules. Specifically, if the coordinate distance between two adjacent species_specific windows on the reference genome of the target species is less than a preset max_gap threshold (preferably 50~200 nt), then the two and their intervening regions are merged into a single target; the merged target is processed uniformly according to its full length when tiling is performed in step S5.

[0042] Preferably, in step S5, the probe length is 80~150 nt, preferably 120 nt; the probe step size is 1 / 5~1 / 3 of the probe length, preferably 1 / 5 of the probe length.

[0043] Preferably, after step S5 and before step S6, a step of classifying host off-target effects for each candidate probe is included: performing a search on the host reference genome for each candidate probe to obtain the nuclear mitochondrial transposon insertion fragment annotation track, protein-coding gene exon annotation, and pseudogene annotation of the host reference genome; classifying each off-target hit into four categories: confirmed, pseudo, borderline, and clean; wherein confirmed probes are removed from the panel, pseudo probes are retained, and borderline probes are entered into the manual review queue. Specifically, a BLAST search is performed on the host reference genome for each candidate probe to obtain the off-target hit set for each probe on the host reference; the nuclear mitochondrial transposon insertion fragment (NUMT) annotation track of the host reference genome is obtained, which contains the coordinates of the annotated mitochondrial transposon insertion segments in the host nuclear genome; protein-coding gene exon annotations and annotated pseudogene annotations of the host reference genome are obtained; each off-target hit is classified as confirmed, pseudo, borderline, or clean; probes classified as confirmed are removed from the candidate oligonucleotide probe sequence set, probes classified as pseudo are retained, and probes classified as borderline are entered into a manual review queue; and the candidate oligonucleotide probe sequence set is updated based on the removal results. A confirmed classification means that the off-target hit falls on a protein-coding gene exon or transcriptional active site, and its identity is not lower than a preset off-target identity threshold, and its hit length is not lower than a preset off-target length threshold. The classification as "pseudo" means that the off-target hit falls within an annotated pseudogene region or within the region covered by the NUMT annotation track. "Borderline" means that the off-target hit falls within a low-complexity region, an unannotated region, or a gray area between "confirmed" and "pseudo." "Clean" means that the probe has no hits on the host reference or the hits are far below a preset threshold. This invention uses NUMT to sense four levels of host off-target classification, avoiding misclassification of mitochondrial-derived probes as high-risk off-targets and protecting the core species identification function of the panel.

[0044] Preferably, the preset miss identity threshold is 70%~95%, preferably 80%; the preset miss length threshold is 40~100 nt, preferably 60 nt.

[0045] Preferably, before performing tiling in step S5, a fusion step of candidate target regions from other sources is included. These candidate target regions include at least one of the following: species-specific marker sequences, species-related antigen gene sequences, high-copy repeat marker sequences, full-length mitochondrial genome tiling sequences, or other species-specific sites in the nuclear genome. For high-copy repeat marker sequences, search hits that meet a preset threshold are retained on each contig of the target species' reference genome. This invention improves the completeness of panel coverage and the accuracy of depth estimation through multi-source structured target fusion and full remapping of high-copy repeat markers.

[0046] Preferably, when the host off-target classification step and the fusion step of candidate target regions from other sources are performed simultaneously: the candidate targets generated by the multi-source fusion step and the species-specific candidate targets generated by the dual-engine determination are merged to form a unified target set before the host off-target classification step is executed. This ensures that multi-source candidate targets and dual-engine candidate targets are treated consistently during the off-target screening process.

[0047] Preferably, before outputting in step S6, the method further includes dividing the candidate oligonucleotide probe sequence set into two subsets: a Core panel and Research add-ons. This division is based on a comprehensive score Starget, which is calculated by weighting the target's source confidence score, host off-target safety score, and threshold robustness score. Targets with an Starget greater than or equal to a preset Core threshold are assigned to the Core panel, while the rest are assigned to Research add-ons. This invention, through hierarchical output of Core and Research add-ons, allows for flexible configuration of panel design results according to clinical priorities.

[0048] Preferably, after generating the candidate oligonucleotide probe sequence set in step S5 and before outputting in step S6, the method further includes filtering probes according to at least one of the following criteria: removing probes with GC content exceeding a preset range; removing probes with the longest continuous homopolymer length exceeding a preset threshold; and removing probes with a low complexity score exceeding a preset threshold.

[0049] Preferably, step S6 further includes outputting a process audit table, which records at least one of the following fields for each candidate window, each target, and each probe: the name of the alignment engine used (engine_used), the judgment result for each closely related species reference (per_ref_verdict), the maximum identity among multiple references (max_identity), the maximum coverage among multiple references (max_coverage), and whether fallback was triggered (fallback_triggered). This invention outputs traceable audit fields for each judgment step, facilitating process auditing and quality control.

[0050] Preferably, the probe panel is used for hybridization enrichment of capture sequencing libraries for at least one of the following sample types: circulating cell-free DNA (cfDNA), formalin-fixed paraffin-embedded (FFPE) sample-derived DNA, and fresh tissue-derived DNA; the capture sequencing library is compatible with at least one of the following next-generation sequencing platforms: Illumina, MGI, Element Biosciences, Ultima Genomics, and Singular Genomics.

[0051] Preferably, the target species and the closely related species are different species that are taxonomically closely related, and the close relationship is defined by at least one of the following criteria: (i) phylogenetic relationship of the same genus or complex; (ii) average nucleotide similarity (ANI) of the whole genome is not lower than a specific threshold of the taxonomic level, of which bacteria are not lower than 95%, fungi are not lower than 80%, and parasites are not lower than 70%; (iii) nuclear genome collinearity block coverage is not lower than 60%; (iv) the most recent common ancestor relationship in a recognized taxonomic database (such as NCBI Taxonomy, GTDB); the target species and the closely related species are preferably selected from different species within the same genus of parasites, fungi, and bacteria.

[0052] Preferably, the dual-engine adaptive cross-species comparison in step S3 adopts at least one of the following computational acceleration methods: the main engine and the fallback engine are executed in parallel on different CPU cores or different computing nodes; the candidate window set is divided into segments and distributed to multiple nodes for parallel processing; the main engine or the fallback engine is implemented using GPU acceleration.

[0053] Preferably, the method further includes processing consecutive N-base pairs (unsequencing gap sequences) present in the reference genomes of the target species and the closely related species according to the following rules, wherein the processing includes at least one of the following:

[0054] (a) Target reference contains N: When the candidate window sequence cut out in step S2 contains N bases, at least one of the following shall be adopted: the candidate window sequence shall be excluded from the output set of step 2; the candidate window sequence shall be divided into two or more sub-windows that do not contain N and sent to step S3 for comparison respectively, and the judgment results obtained shall be merged according to the conservative rule that any sub-window triggers shared, i.e. the whole is shared.

[0055] (b) N-related scenarios in closely related references: When the reference genome of the closely related species is a scaffold or contig-level unchromosome-assembled genome, and the alignment result of the reference by the main engine or fallback engine in step S3 spans its gap sequence, the N base positions covered by the gap sequence are excluded when calculating the optimal local alignment segment length (region_len) and the number of matching bases, so as to avoid falsely low pident due to N bases being miscalculated as mismatch or falsely high region_len due to N bases being miscalculated as match.

[0056] This invention also provides a capture sequencing probe panel design system for the genomes of closely related species, comprising:

[0057] The reference genome management module is used to obtain the reference genome sequence of the target species and the reference genome sequence of closely related species.

[0058] The candidate window generation module is used to extract candidate window sequences with a length of no more than 200 nucleotides from the reference genome of the target species at a fixed step size;

[0059] a dual-engine adaptive cross-species alignment module, configured to perform alignment with a main engine based on a seed-extension exact alignment algorithm, obtain a first optimal local alignment segment satisfying a preset e-value threshold, calculate the fields of percent identity and query coverage, and record them as (pident_p, qcovs_p); wherein the percent identity field is defined as the ratio of the number of matched bases in the first optimal local alignment segment to the length of the first optimal local alignment segment, and the query coverage field is defined as the ratio of the length of the first optimal local alignment segment to the length of a candidate window sequence; when any preset trigger condition is satisfied, a fallback engine based on the preset assembly-level low difference degree is triggered to perform re-alignment, a second optimal local alignment segment is selected, and percent identity and query coverage are calculated according to the same field definition rule as that of the main engine, and recorded as (pident_f, qcovs_f); if the first optimal local alignment segment is obtained and any trigger condition is not satisfied, (pident_p, qcovs_p) is used to enter a fusion determination module; otherwise, (pident_f, qcovs_f) is used to enter the fusion determination module;

[0060] a fusion determination module, configured to record (pident_p, qcovs_p) or (pident_f, qcovs_f) output by the dual-engine adaptive cross-species alignment module as (identity, coverage), and perform the following determination on each candidate window based thereon: if the determination for any closely-related species reference genome satisfies identity≥T_shared and coverage≥T_cov, the fusion determination result is shared; otherwise, if the determination for any closely-related species reference genome satisfies T_cross≤identity<T_shared and coverage≥T_cov, the fusion determination result is cross_reactive; in other cases, the fusion determination result is species_specific; wherein T_shared, T_cross and T_cov are preset thresholds;

[0061] a probe tiling module, configured to collect candidate windows determined as species_specific by the fusion determination module as targets, and perform tiling to generate a candidate oligonucleotide probe sequence set;

[0062] a probe set output module, configured to output the candidate oligonucleotide probe sequence set as a species-specific capture probe panel.

[0063] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements any one of the methods described herein.

[0064] The present invention also provides a computing device, including a processor and a memory, the memory storing a computer program, wherein the processor executes the computer program to implement any one of the methods described above.

[0065] The beneficial effects of this invention are:

[0066] (1) This invention strictly defines the calculation rules for the decision fields pident and qcovs, calculating only based on the optimal local alignment segment itself, thus excluding the alignment block length as the basis for coverage estimation from the source. This avoids the systematic false positives caused by anchor point extension in short query scenarios by alignment tools based on short read presets. In short query close species alignment scenarios with no more than 200 nucleotides, the false positive rate of species-specific window determination in this invention is reduced from about 45% to 65% to less than 5% compared to the traditional single-engine scheme based on short read presets.

[0067] (2) This invention provides fallback engine triggered by conditions to provide fallback judgment in edge scenarios such as main engine NO-HIT, low bit score, short HSP, and high e-value, avoiding the judgment loss of single engine solution in these scenarios and restoring about 8% to 12% of the decision window.

[0068] (3) This invention retains three classification levels: shared, cross_reactive, and species_specific, and provides quantifiable intermediate criteria for 70% to 90% identity edge cases that are common among closely related species. Compared with the binary judgment scheme, it has higher discrimination accuracy. Attached Figure Description

[0069] Figure 1 This is a flowchart of the dual-engine adaptive design method provided in the embodiments of the present invention;

[0070] Figure 2 This diagram illustrates the failure mechanism of existing alignment tools based on short read presets in scenarios with short queries of no more than 200 nucleotides.

[0071] Figure 3 This is a flowchart of the dual-engine adaptive decision-making process provided in an embodiment of the present invention;

[0072] Figure 4 This is a comparison chart of the species-specific ratios between the dual-engine determination method of Embodiment 1 of the present invention and existing methods;

[0073] Figure 5 This is a comparison diagram of two alignment methods under the same 200 nucleotide query window in Embodiment 1 of the present invention. Detailed Implementation

[0074] The following specific embodiments illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0075] refer to Figure 1 The dual-engine adaptive design method for short query of capture sequencing probe panels for closely related species genomes provided in this invention mainly includes the following steps: obtaining the target species and reference genomes of closely related species; cutting out candidate window sequences no larger than 200 nt; performing dual-engine adaptive cross-species alignment for each candidate window sequence, where the main engine selects the first optimal local alignment segment that meets the preset e-value threshold based on the seed-extended precise alignment algorithm and calculates the pident and qcovs fields; when a preset trigger condition is met, the condition-triggered fallback engine selects the second optimal local alignment segment and calculates according to the same field definition rules; performing single-reference three-level judgment and multi-reference fusion judgment based on identity and coverage; collecting species_specific windows as the target set and performing tiling to generate a candidate probe sequence set. The following provides a detailed description of each step.

[0076] Step S1: Reference Genome Acquisition and Preprocessing

[0077] Obtain the reference genome sequence of the target species and one or more reference genome sequences of closely related species. The reference genome can be a complete assembled sequence in FASTA format, or an assembled sequence at the contig or scaffold level. For the target species, its genome size is typically greater than or equal to 1 Mb. The number of closely related species is preferably 2 to 4.

[0078] Step S2: Cut out the candidate window sequence

[0079] Candidate window sequences are cut from the reference genome of the target species at a fixed step size and window size. The length of each candidate window sequence is no more than 200 nucleotides, preferably 120 to 200 nucleotides. The step size is preferably between one-half and one-fifth of the window size; for example, when the window size is 200 nucleotides, the step size can be 50 to 100 nucleotides.

[0080] The basic process for cutting out candidate windows is as follows: starting from the 5' end of the reference sequence, subsequences of length windowsize are successively cut as candidate windows; after each cut, the starting position is moved backward by step length; if the remaining length at the end of the sequence is less than window size, the remaining segment is used as the last candidate window.

[0081] Step S3: Dual-engine adaptive cross-species alignment

[0082] refer to Figure 3 Step S3 is the core step of this invention. For each candidate window and for each closely related species' reference genome, a dual-engine adaptive cross-species alignment is performed. This step includes sub-steps S3a (main engine alignment), S3b (fallback engine conditional triggering), and S3c (result selection and output). The overall decision-making process is as follows: Figure 3 As shown.

[0083] Sub-step S3a: Main engine alignment. The candidate window is aligned with the reference genome of closely related species using a seed-and-extend precise alignment algorithm. The main engine is selected from at least one of blastn, megablast, blastn-short, BLAT, LAST, USEARCH, edlib, or parasail, with blastn being preferred. The blastn parameters are set as follows: word_size between 7 and 15 (11 recommended for short queries), e-value between 1e-3 and 1e-8 (1e-5 recommended), dust filtering is disabled, and gapopen and gapextend use their default values. After main engine alignment, the best high-scoring pair (best HSP) that meets the preset e-value threshold is selected. For the best HSP, two key fields are calculated: the percent identity field (pident_p) and the query coverage field (qcovs_p).

[0084] The calculation rule for pident_p is: the number of matching bases (identities) in the best HSP divided by the length of the best HSP (alignment_length), that is:

[0085] pident_p = hsp.identities / hsp.alignment_length * 100

[0086] The calculation rule for qcovs_p is: the best HSP length (alignment_length) divided by the candidate window sequence length (query_length), that is:

[0087] qcovs_p = hsp.alignment_length / query_length * 100

[0088] The two fields mentioned above are calculated solely based on the best HSP and not on any extended alignment block. This is one of the core differences between this invention and existing technologies: existing methods often calculate coverage based on the entire alignment block, leading to an artificial overestimation of coverage in short query scenarios; this invention strictly limits the use of only the length of the best HSP, eliminating this systematic bias from the source.

[0089] It is understood that the term "optimal local alignment segment" is a general concept across tools. Table 1 below shows the specific correspondence between the optimal local alignment segments in the original output formats of several preferred main engine tools of this invention.

[0090] Table 1. Correspondence rules of the optimal local alignment segment in each main engine tool

[0091]

[0092] In the above correspondence rules, the number of matching bases uniformly refers to the number of bases that are strictly identical to the target reference (excluding mismatch, gap, and N bases); the region length (region_len) refers to the coverage length of the optimal local alignment region on the query sequence (including mismatch but excluding query-end gap). Through this unified definition, the two decision fields, pident=number of matching bases / region_len and qcovs=region_len / query_len, maintain semantic consistency across all main engine tools.

[0093] Sub-step S3b: Fallback engine conditional triggering. When at least one of the following conditions is met, the conditional triggering invokes the fallback engine based on the assembly-level low-difference preset for re-alignment:

[0094] (i) The main engine did not output any optimal local alignment segment (NO-HIT);

[0095] (ii) The maximum bit score of the optimal local comparison segment output by the main engine is lower than the preset lower limit of bit score (30 to 50 recommended).

[0096] (iii) The maximum segment length of the optimal local alignment segment output by the main engine is lower than the preset length limit (20 to 30 nucleotides recommended).

[0097] (iv) The minimum e-value of the optimal local comparison segment output by the main engine is greater than the preset e-value threshold.

[0098] The fallback engine is selected from at least one of minimap2 -cx asm5, -cx asm10, -cx asm20, LAST, or nucmer, preferably minimap2 -cx asm10. The fallback engine is suitable for scenarios where the main engine cannot provide a valid HSP due to excessively high differences between the short query window and the reference genome of closely related species, or where the HSP provided by the main engine is at the edge of confidence. After the fallback engine outputs, the best HSP is extracted and pident_f and qcovs_f are obtained according to the same calculation rules as sub-step S3a. In some implementations, the two engines can be executed serially (whether to call fallback depends on the trigger condition after the main engine completes) or in parallel (the main engine and fallback engine run simultaneously; if the main engine returns a valid result, the result from the main engine is used; otherwise, the result from the fallback engine is used). Parallel mode is suitable for scenarios with high concurrent computing resources to shorten end-to-end time.

[0099] Sub-step S3c: Result selection and output. If sub-step S3a gives a valid result and does not meet any of the triggering conditions described in sub-step S3b, then use pident_p and qcovs_p output by sub-step S3a to proceed to step S4; otherwise, use pident_f and qcovs_f output by the fallback engine; if neither sub-steps S3a nor S3b gives a valid result, then the relationship between the window and the reference is determined to be no significant homology, and the identity and coverage fields of the window-reference pair are recorded as missing values.

[0100] Step S4: Single-reference three-level determination and multi-reference fusion determination

[0101] Based on the identity (abbreviated as id) and coverage (abbreviated as cov) fields output in step S3, a single-reference three-level judgment is performed on each candidate window, followed by a multi-reference fusion judgment.

[0102] The single reference three-level determination rules are as follows: For a candidate window with id and cov values ​​for a reference genome of a closely related species: if id is greater than or equal to T_shared and cov is greater than or equal to T_cov, then the single reference is determined as shared; if id is greater than or equal to T_cross and less than T_shared and cov is greater than or equal to T_cov, then the single reference is determined as cross_reactive; otherwise, the single reference is determined as species_specific.

[0103] The multi-reference fusion determination rules are as follows: A candidate window typically needs to be aligned with reference genomes of multiple closely related species. Fusion determination employs a trigger-and-upgrade strategy: if the window is determined as shared for any single reference of a closely related species, then the fusion determination is shared; otherwise, if the window is determined as cross-reactive for any single reference, then the fusion determination is cross-reactive; the fusion determination is only specific if the single reference determination for all references is species-specific.

[0104] The preferred preset thresholds are: T_shared 85% to 95% (90% recommended), T_cross 65% to 80% (70% recommended), and T_cov 70% to 90% (80% recommended). These thresholds can be adjusted based on the actual evolutionary distance between the target species and closely related species.

[0105] Step S5: Target Collection and Tiling Probe Generation

[0106] Collect all candidate windows whose fusion is determined to be species_specific as a target set. Perform tiling on each target in the target set according to the preset probe length (probe_length) and probe step (probe_step) to generate a set of candidate oligonucleotide probe sequences.

[0107] The basic procedure for tiling is as follows: For each target sequence, starting from the 5' end, sequentially extract subsequences of length probe_length as candidate probes; after each extraction, shift the starting position forward by probe_step; if the remaining length at the end of the sequence is less than probe_length, then that end fragment is used as the last candidate probe. The probe length is preferably 80 to 150 nucleotides (120 nucleotides recommended), and the probe step size is preferably between one-fifth and one-third of the probe length (24 nucleotides recommended, corresponding to one-fifth of the 120 nucleotide probe length).

[0108] Before tiling, candidate target regions from other sources can be fused with the target set obtained in step S4. Other candidate target regions include: species-specific marker sequences reported in the literature, species-related antigen gene sequences, high-copy repeat marker sequences, full-length mitochondrial genome tiling sequences, or other species-specific sites in the nuclear genome. Multi-source fusion occurs before tiling; candidate targets from all sources undergo a unified downstream quality control process before being uniformly tiled.

[0109] Step S6: Panel Output

[0110] The output consists of a set of candidate oligonucleotide probe sequences as a species-specific capture probe panel. The output preferably includes at least one of the following: a list of nucleotide sequences of candidate probes, target source information for each probe, the probe's coordinates on the target genome, and a process audit table. The process audit table records the alignment engine name used for each candidate window, each target, and each probe; the single-reference three-level judgment result for each closely related species reference; the maximum identity and maximum coverage among multiple references; and whether fallback was triggered, enabling full traceability of the process.

[0111] Before output, one or more of the following quality controls can be performed on the candidate probe sequence set:

[0112] (a) Host off-target classification: For each candidate probe, a search is performed on the host reference genome to obtain nuclear mitochondrial transposon insertion (NUMT) annotation tracks, protein-coding gene exon annotations, and pseudogene annotations, which are classified into four categories: confirmed, pseudo, borderline, and clean. Confirmed probes are removed from the panel, pseudo probes are retained, and borderline probes are entered into the manual review queue.

[0113] (b) Threshold sensitivity scan: Species-specific determination is performed on the pre-target set at a set of candidate thresholds, and a target is included in the final target set only if the determination is consistent at all thresholds.

[0114] (c) Core / Research layering: Based on the comprehensive score Starget (which is calculated by weighting source confidence score, host off-target safety score and threshold robustness score), probes are divided into two subsets: Core panel and Research add-on.

[0115] (d) Probe filtering: Remove probes with GC content exceeding the preset range (preferably 25% to 65%), longest continuous homopolymer length exceeding the preset threshold (preferably greater than or equal to 8 nucleotides), or low complexity score exceeding the preset threshold.

[0116] Example 1: Species-specific window screening based on reference sequences of closely related eukaryotic pathogens

[0117] Input: Reference sequences of three groups of closely related eukaryotic pathogens, denoted as batch A, batch B, and batch C, respectively. The length of a single sequence in each batch is between 10 and 30 kb. The candidate window generation parameters are: window size equal to 200 nucleotides and step equal to 100 nucleotides, with approximately 100 to 200 candidate windows per batch.

[0118] Control group 1 adopted a single-engine approach based on short-read presets: the minimap2 short-read preset was used to perform alignments against other references within the batch for each candidate window. The alignment block-level identity and the ratio of alignment block to query length from the PAF output were parsed as query coverage fields, and the criteria were determined based on thresholds of T_shared equaling 90%, T_cross equaling 70%, and T_cov equaling 80%. The results showed that the proportion of species_specific windows in all three batches was approximately 45% to 65%.

[0119] Control group 2 employed a BLASTn single-engine plus three-threshold scheme: BLASTn was used to perform alignments for each candidate window against other references within the batch, judging based on the best HSP's pident and qcovs fields and the same three thresholds, without invoking any fallback engine. The results showed that the species_specific proportions for windows where the main engine could provide results were close to the results of this invention. However, BLASTn returned NO-HIT for approximately 8% to 12% of windows per batch. These windows could only be classified as having no significant homology or excluded in the single-engine scheme, affecting the completeness of downstream target collection.

[0120] This invention employs a dual-engine adaptive scheme: blastn is used as the main engine with parameters word_size equal to 11 and evalue equal to 1e-5. For each candidate window, alignment is performed against other references within the batch. When the main engine meets any triggering condition on the query-reference pair, minimap2 -cx asm10 is automatically called as the fallback engine. A threshold is used for judgment. The results show that the proportion of species_specific windows in all three batches is less than 5%, the shared proportion is approximately 25% to 35%, and the cross_reactive proportion is approximately 60% to 70%. (Reference) Figure 4 Compared to a single short-read preset scheme, the present invention reduces false positives by 40 to 60 percentage points; compared to a single BLASTn scheme, it recovers an additional 8% to 12% of the decisionable window in NO-HIT edge scenarios.

[0121] Sequence comparison evidence for determining accuracy: Taking a typical 200-nucleotide query window from Example 1, the alignment block output by the alignment tool based on short read presets and the best HSP output by blastn were extracted respectively. See details for comparison. Figure 5 The alignment tool based on short-read presets outputs an alignment block length of approximately 160 nucleotides, with an apparent identity of approximately 88% to 92%, triggering a shared determination; while blastn reports a best HSP length of approximately 28 nucleotides, with a pidentity of approximately 92% to 95%, correctly identifying it as species-specific. The difference between the two stems from different assumptions made by the algorithms regarding the extension behavior of local anchor points.

[0122] Fallback triggering tests: In each batch of approximately 100 to 200 candidate windows, the blastn main engine triggered fallback for at least one reference window in about 8% to 12% of the windows, of which about 6% to 9% were no-hit, and the rest were low-bit-score or short-HSP edge scenarios. The fallback engine provided quantification criteria for these windows, with about 60% falling in the cross_reactive category, about 30% in the species_specific category, and about 10% in the shared category.

[0123] Example 2: NUMT-based host off-target classification

[0124] Input: The target set (approximately 2000 targets) obtained from the dual-engine determination in Example 1, after tiling, forms a candidate probe set (approximately 70,000 120 nt probes), GRCh38 human reference, UCSC nuMtSeq track BED file, and GENCODE annotation GTF.

[0125] Steps: Perform a host reference genome search for each candidate probe; for each search hit, check if its coordinates fall within the nuMtSeq track segment, GENCODE protein-coding exon, annotated pseudogene, or other segment; classify according to four criteria, using an identity threshold of 80% and a length threshold of 60 nucleotides. Probes classified as "confirmed" are removed from the candidate oligonucleotide probe sequence set, probes classified as "pseudo" are retained, and probes classified as "borderline" are added to the manual review queue; and the candidate oligonucleotide probe sequence set is updated based on the removal results. "Confirmed" means the off-target hit falls on a protein-coding gene exon or transcriptional active site, and its identity is not lower than the preset off-target identity threshold, and its hit length is not lower than the preset off-target length threshold. "Pseudo" means the off-target hit falls within an annotated pseudogene segment or within the segment covered by the NUMT annotation track. "Borderline" means the off-target hit falls within a low-complexity segment, an unannotated segment, or a gray area between "confirmed" and "pseudo". The term "clean" refers to the probe having no hits on the host reference or having hits far below a preset threshold.

[0126] Results: The four-level classification distribution was approximately 0 confirmed, 62 pseudo, 14 borderline, and 2166 clean.

[0127] Control experiment (without NUMT sensing): If the same data were roughly categorized by simple hit count, all approximately 180 probes from mtDNA sources would be misclassified as risk and removed, compromising the core of the panel's species identification. However, the NUMT sensing classification of this invention correctly retains these probes as pseudo files and includes them in the final panel.

[0128] Example 3: Threshold Sensitivity Scan

[0129] Input: Approximately 2000 candidate targets obtained based on the dual-engine determination in Example 1.

[0130] Operation: Determine the species-specific verdict for each target at threshold sets of 10%, 20%, 25%, 30%, 40%, and 50%.

[0131] Output: A two-dimensional sensitivity matrix of target and threshold, approximately 2000 multiplied by 6, with a total of 12000 verdict values.

[0132] Enhanced decision-making: Targets are included in the species-specific target whitelist only if their verdicts are consistent across six thresholds. In this example, approximately 785 targets passed the consistency test, while the remaining targets were marked as threshold-sensitive.

[0133] Example 4: Full remapping of high-copy repeat marker

[0134] Input: A repeat marker sequence, the target species' reference genome contains approximately 1200 contigs.

[0135] Compared to the industry default approach, blastn retains the first 5 hits for each query by default, but in reality only captures 5 copy locations.

[0136] The embodiments of the present invention employ full contig search to capture all approximately 600 copy sites of the marker in the genome, and generate corresponding tiling probes for each copy site to ensure accurate actual panel coverage depth.

[0137] Example 5: Complete panel design and end-to-end operation

[0138] Complete end-to-end operational parameters for the experimental process: Input of one target species reference genome, two to four closely related species reference genomes, and multiple types of multi-source markers. Key engineering parameters are: probe length 120 nucleotides, probe step size 24 nucleotides, and panel limit 4 Mb.

[0139] The output panel has approximately 2200 targets, a total length of approximately 2Mb, and approximately 72000 probes. After separating the Core and Research panels, approximately 350 targets are assigned to the Core panel, and approximately 1900 targets are assigned to the Research add-on.

[0140] The runtime is approximately 4 to 6 hours, but can be shortened to approximately 1.5 hours by enabling parallel execution of the main engine and fallback engine, as well as multi-node parallel execution with candidate window sharding.

[0141] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope of the present invention.

Claims

1. A short query design method for capture sequencing probe panels, characterized in that, Comprising: S1. Obtaining a reference genome sequence of a target species and reference genome sequences of closely related species; S2. Cutting out candidate window sequences with a length of no more than 200 nucleotides from the reference genome of the target species at a fixed step size; S3. Performing dual-engine alignment on the candidate windows and the reference genomes of closely related species, comprising: S3a. Performing alignment with a main engine based on a seed-extension exact alignment algorithm, obtaining a first optimal local alignment segment that satisfies a preset e-value threshold, and calculating the fields percent identity and query coverage, recorded as (pident_p, qcovs_p); wherein the percent identity field is defined as the ratio of the number of matched bases in the optimal local alignment segment to the length of the optimal local alignment segment, and the query coverage field is defined as the ratio of the length of the optimal local alignment segment to the length of the candidate window sequence; S3b. When any preset trigger condition is satisfied, triggering a fallback engine based on the assembly-level low difference preset to perform realignment, selecting a second optimal local alignment segment, calculating percent identity and query coverage according to the same field definition rule as the main engine, and recording the calculation result as (pident_f, qcovs_f); S3c. If the first optimal local alignment segment is obtained in sub-step S3a and any trigger condition in sub-step S3b is not satisfied, step S4 is performed using (pident_p, qcovs_p); otherwise, step S4 is performed using (pident_f, qcovs_f) of sub-step S3b; S4. Recording (pident_p, qcovs_p) or (pident_f, qcovs_f) output in step S3 as (identity, coverage), and making the following determination for each candidate window based thereon: if the determination for any reference genome of a closely related species satisfies identity≥T_shared and coverage≥T_cov, the fusion determination is shared; otherwise, if the determination for any reference genome of a closely related species satisfies T_cross≤identity<T_shared and coverage≥T_cov, the fusion determination is cross_reactive; in other cases, the fusion determination is species_specific; wherein T_shared, T_cross and T_cov are preset thresholds; S5. Collecting candidate windows whose fusion determination is species_specific in step S4 as targets, and performing tiling to generate a candidate oligonucleotide probe sequence set; S6. Outputting the candidate oligonucleotide probe sequence set as a species-specific capture probe panel.

2. The method as described in claim 1, characterized in that, In step S3, The selection rule for the first optimal local alignment segment is as follows: among all local alignment results output by the main engine on the candidate window and the reference genome pair of closely related species, the local alignment segment with the maximum bit score is selected as the first optimal local alignment segment; and / or, The selection rule for the second optimal local alignment segment is as follows: among all local alignment results output by the fallback engine on the candidate window and the reference genome pair of closely related species, the local alignment segment with the maximum value of the actual matching segment length is selected as the second optimal local alignment segment; and / or, The percent identity and query coverage fields are calculated based on the actual matching segments within the selected first or second best local alignment segment, rather than using the total length of the alignment block as the basis for coverage estimation.

3. The method as described in claim 1, characterized in that, In sub-step S3b, the preset triggering condition includes at least one of the following: (i) The main engine did not output the first optimal local alignment segment; (ii) The maximum bit score of the first optimal local alignment segment is lower than the preset lower limit of the bit score; (iii) The maximum segment length of the first optimal local comparison segment is lower than the preset length lower limit; (iv) The minimum e-value of the first optimal local alignment segment is greater than the preset e-value threshold; The preset lower limit for bit score is 30-50, the preset lower limit for length is 20-30 nucleotides, and the preset e-value threshold is 10. -3 ~10 -8 .

4. The method as described in claim 1, characterized in that, If neither sub-steps S3a nor S3b yields the first optimal local alignment segment nor the second optimal local alignment segment, then the candidate window is determined to have no significant homology with the reference genome of the closely related species.

5. The method as described in claim 1, characterized in that, In sub-step S3a, the main engine is selected from at least one of blastn, megablast, blastn-short, BLAT, LAST, USEARCH, edlib, and parasail; and / or, In sub-step S3b, the fallback engine is selected from at least one of minimap2 -cx asm5, -cx asm10, -cx asm20, LAST, and nucmer; and / or, In step S4, T_shared is 85%~95%, T_cross is 65%~80%, and T_cov is 70%~90%.

6. The method as described in claim 1, characterized in that, After step S5 and before step S6, the method further includes a step of classifying host off-target effects for each candidate probe: performing a search on the host reference genome for each candidate probe to obtain the nuclear mitochondrial transposon insertion fragment annotation track, protein-coding gene exon annotation, and pseudogene annotation of the host reference genome; classifying each off-target hit into four categories: confirmed, pseudo, borderline, and clean; wherein confirmed probes are removed from the panel, pseudo probes are retained, and borderline probes are entered into the manual review queue.

7. The method as described in claim 1, characterized in that, After step S4 and before step S5, a threshold sensitivity scanning step is also included: for the pretarget set composed of candidate windows that are determined to be species_specific by fusion, each pretarget is determined to have species specificity under a set of preset candidate thresholds, and a two-dimensional sensitivity matrix of target and threshold is output. A target is included in step S5 only if the results of a pre-target are consistent across all candidate thresholds.

8. The method as described in claim 1, characterized in that, Before tiling in step S5, there is also a fusion step of candidate target regions from other sources. The candidate target regions from other sources include at least one of species-specific marker sequences, species-related antigen gene sequences, high-copy repeat marker sequences, full-length mitochondrial genome tiling sequences, or other species-specific sites in the nuclear genome. For high-copy repeat marker sequences, the search hits that meet the preset threshold on each contig of the target species reference genome are retained.

9. The method as described in claim 1, characterized in that, Before the output of step S6, the method further includes dividing the candidate oligonucleotide probe sequence set into two subsets: Core panel and Research add-on. The division is based on a comprehensive score Starget, which is calculated by weighting the source confidence score, host off-target safety score, and threshold robustness score of the target. Targets with Starget greater than or equal to a preset Core threshold are assigned to the Core panel, and the rest are assigned to Research add-on.

10. The method as described in claim 1, characterized in that, After generating the candidate oligonucleotide probe sequence set in step S5 and before outputting in step S6, the method further includes filtering probes according to at least one of the following criteria: removing probes with GC content exceeding a preset range; removing probes with the longest continuous homopolymer length exceeding a preset threshold; and removing probes with a low complexity score exceeding a preset threshold.

Images