Design method of clinical metagenome Spike-in sequence

Abstract

The invention discloses a design method of a clinical metagenome Spike-in sequence, and belongs to the field of metagenome detection. The method comprises the following steps: by taking an arabidopsis thaliana genome as a sequence source, firstly segmenting the arabidopsis thaliana genome into 1000bp candidate fragments, and then splitting the candidate fragments into 50bp k-mer; the method comprises the following steps: strictly screening a sequence without microbe and human homologous k-mer through Kraken2 conservative parameters (the minimum k-mer matching length is 31bp, the matching fraction is 0.95 and the like), then carrying out specific screening of'non-repeated k-mer 'or'minimum editing distance is greater than or equal to 10', and finally carrying out homology and physicochemical property verification to obtain a final Spike-in sequence. The problems that an existing Spike-in sequence is prone to being homologous with a microorganism sequence or a human source sequence, and k-mer distinguishability among multiple sequences is insufficient are solved, the designed Spike-in sequence is high in safety and distinguishability and adapts to a clinical sequencing process, synchronous quality control and accurate absolute quantification of multiple internal standards can be achieved, the accuracy of clinical pathogenic microorganism detection is improved, and the method is suitable for popularization and application. The method is suitable for clinical diagnosis scenes of difficult infection and critical infection.

Description

Technical Field

[0001] This invention relates to metagenomic detection, and more particularly to a method for designing clinical metagenomic Spike-in sequences. Background Technology

[0002] Metagenomic high-throughput sequencing (mNGS) technology, with its advantages of no culture required and wide coverage, has become a core tool for the diagnosis of difficult and critical infections in clinical practice, enabling rapid identification of various pathogens such as bacteria, fungi, and viruses in samples. However, the clinical metagenomic sequencing process involves multiple steps, including nucleic acid extraction, PCR amplification, library construction, sequencing, and data analysis. It is susceptible to factors such as technical bias (e.g., differences in extraction efficiency, amplification bias) and sample complexity (e.g., high proportion of host nucleic acid, low pathogen abundance), leading to quantitative deviations or false positives / false negatives in sequencing results.

[0003] Spike-in technology (exogenous internal standard technology) is a key means to solve the above problems—by introducing exogenous standards with known sequences and copy numbers into clinical samples, it participates in the sequencing process throughout the entire process, realizing full-process quality control, cross-batch data correction, and absolute quantification of pathogens. However, there are two major technical bottlenecks in the design of Spike-in sequences: (1) Risk of microbial homology interference: If the Spike-in sequence is homologous to the microbial sequence that may exist in the clinical sample, it will be misjudged as a microorganism by species annotation software (such as Kraken2), interfering with the detection results of the real pathogen. (2) Insufficient sequence distinguishability: In clinical testing, multiple Spike-in sequences (such as internal standards with different concentration gradients) are often used at the same time. If the short fragments (k-mers) of each sequence are repetitive or highly similar, it will be impossible to accurately trace the original Spike-in sequence corresponding to each k-mer after sequencing, thus affecting the accuracy of internal standard counting and the reliability of multi-internal standard simultaneous analysis.

[0004] Existing Spike-in sequence designs mostly use artificially synthesized random sequences or select rare microbial genome fragments, which have the following drawbacks: (1) The safety and stability of random sequences are difficult to guarantee, and homology with unknown microbial sequences cannot be completely ruled out; (2) Rare microbial genome fragments still have the potential risk of microbial homologous sequences, and the sources are limited and difficult to prepare on a large scale; (3) There is a lack of standardized homology screening and sequence-specific screening procedures, which can not accurately filter microbial homologous sequences, nor consider the k-mer distinguishability between multiple sequences, resulting in internal standards interfering with the identification of target pathogens or multiple internal standards being confused in clinical testing.

[0005] Arabidopsis thaliana, as a model organism, has a fully resolved genome sequence and exhibits extremely low homology with the genomes of common clinical pathogens. Furthermore, it poses no risk of infecting humans, making it an ideal source of Spike-in sequences. However, directly extracting fragments from the Arabidopsis thaliana genome as Spike-in sequences may still lead to misinterpretation due to accidental homology between short local fragments and microbial sequences. Additionally, multiple sequences may exhibit k-mer repetitions or high similarity, affecting the effectiveness of using multiple internal standards. Summary of the Invention

[0006] Purpose of the invention: The purpose of this invention is to provide a method for designing clinical metagenomic Spike-in sequences, aiming to solve the dual problems of existing Spike-in sequences being easily homologous to microbial or human sequences and insufficient k-mer distinguishability among multiple sequences.

[0007] Technical solution: A method for designing clinical metagenomic Spike-in sequences, comprising the following steps:

[0008] S1. Construction of candidate sequence library

[0009] The Arabidopsis genome (preferably TAIR10 or a newer version of the reference genome) was selected as the sequence source. Sequence extraction tools (such as Bedtools, seqkit, or a custom Python script) were used to divide the Arabidopsis whole genome sequence into candidate fragments of uniform length. The fragment length was set to 1000 bp. During the extraction process, it was ensured that there was no overlap or omission, and the entire Arabidopsis genome was covered to form a candidate sequence library.

[0010] S2, k-mer splitting

[0011] For each 1000bp candidate sequence in the candidate sequence library, a k-mer splitting algorithm (such as a custom Python script or unikmer software) is used to split it into consecutive short k-mer fragments. The k-mer length is set to 50bp, and the splitting rule is: adjacent k-mers overlap by 49bp (i.e., splitting by sliding one base at a time), ensuring that each 1000bp candidate sequence can be split into 951 (1000-50+1) 50bp k-mers, forming the k-mer set corresponding to each candidate sequence.

[0012] S3. Strict screening based on microbial and human homology (first screening)

[0013] Kraken2 software was used to annotate all k-mers for each candidate sequence by species classification. Conservative parameters were set to maximize the filtering of candidate sequences containing potential microbial or human homologous k-mers. The core purpose of this step is to eliminate the possibility of Spike-in sequences in subsequent metagenomic sequencing data being misidentified as human or pathogenic microorganism sequences, thereby completely avoiding their interference with the normal detection of human and pathogenic microorganism sequences in clinical metagenomic samples and ensuring the accuracy and specificity of subsequent pathogen detection results.

[0014] The specific process is as follows:

[0015] (1) Kraken2 database configuration: Select the whole genome database and human genome database that includes all common clinical pathogens such as bacteria, fungi, viruses, mycoplasma, chlamydia, and rickettsia in NCBI RefSeq, and the database is updated regularly to cover new pathogens;

[0016] (2) Kraken2 conservative parameter settings: set the minimum k-mer matching length to 31bp, the minimum matching score to 0.95, turn off fast mode (disabled by the --fast parameter), and enable strict classification mode (the --strict parameter) to ensure that only k-mers with high confidence matching microorganisms or human origins will be labeled;

[0017] (3) Screening criteria: For each candidate sequence, if any one of the 951 50bp k-mers is assigned to a microbial species (i.e., the classification result is bacteria, fungi, viruses, etc.) or human origin by Kraken2, then the candidate sequence is judged as a "potential interference sequence" and filtered; only candidate sequences in which all k-mers are not assigned to a microbial species or human origin by Kraken2 are retained to form the first set of qualified sequences.

[0018] S4. Sequence-specific enhanced screening (secondary screening)

[0019] The final Spike-in sequences are selected from the first set of qualified sequences. The core requirement is that there are no duplicates or the edit distance between each pair of 50bp k-mers in each sequence is as large as possible, ensuring that the original Spike-in sequence corresponding to each k-mer can be accurately traced after sequencing. The core purpose of this step is to ensure that after metagenomic sequencing, all successfully detected Spike-in-related sequences can be accurately located to their corresponding original Spike-in sequences. This design provides a key sequence anchoring basis for subsequent source tracing analysis of samples, effectively ensuring the uniqueness and accuracy of the source tracing results, and laying a solid foundation for the quality control and result verification of the entire detection system.

[0020] The specific screening process is as follows:

[0021] (1) Construct a global k-mer set: Extract the 50bp k-mer of all sequences in the first-stage qualified sequence set to form a global k-mer set, and record the sequence ID corresponding to each k-mer;

[0022] (2) Selection of screening strategies (choose one or use in combination): Strategy A (no duplicate k-mer screening): Traverse the global k-mer set and remove sequences containing duplicate k-mers (i.e., the same k-mer corresponds to multiple sequence IDs); retain sequences where all k-mers are unique to the sequence and are not shared by any other sequences, forming a subset of no duplicate k-mer sequences; Strategy B (high edit distance screening): Perform pairwise comparisons on the qualified sequence set of the first screening, calculate the edit distance (preferably Hamming distance) between all k-mers in the k-mer set corresponding to any two sequences, and screen sequences that satisfy "among all pairwise k-mer pairs of any two sequences, the minimum edit distance is ≥10"; if multiple sequences need to be screened, use a greedy algorithm to select them step by step: first randomly select 1 sequence as the benchmark, and then select sequences whose k-mer edit distances to the selected sequences are all ≥10, until the required number of sequences is reached;

[0023] (3) Determining the number of sequences: Based on clinical testing needs (such as single internal standard or multi-gradient internal standard), select several sequences (preferably 20-30) from the above subset or screening results to form a second set of qualified sequences.

[0024] S5 and Spike-in sequence verification

[0025] The qualified sequence set from the second round of screening is then validated again to ensure compatibility with clinical metagenomic sequencing workflows.

[0026] (1) Secondary verification of microbial homology: The homology of the sequence with the microbial sequence in the NCBI NR / NT database was compared using the BLASTn tool to ensure that the homology between the sequence and the microbial sequence was less than 70%;

[0027] (2) Human sequence homology verification: BLASTn alignment with the human genome (GRCh38 or later version) is performed to ensure that the homology is less than 50% and to avoid cross-reaction with the host sequence;

[0028] (3) Physicochemical property verification: Analyze the GC content of the sequence and retain the sequence with a GC content between 40% and 60% to ensure compatibility with PCR amplification and sequencing procedures;

[0029] (4) Distinguishing verification: For the finally selected sequences, the k-mer repetition rate (preferably ≤0%) or edit distance (preferably ≥25) between each pair is calculated again to verify that the distinguishability meets the standard.

[0030] Beneficial effects:

[0031] (1) Completely avoid interference from microorganisms and human homology: Through the dual guarantee of "Arabidopsis thaliana genome source + Kraken2 conserved parameter screening", microorganisms and human homology sequences are comprehensively screened at the short fragment (50bp k-mer) level, ensuring that the final Spike-in sequence has no k-mers that can be identified as microorganisms or humans, thus completely avoiding interference from internal standards in the detection and annotation of real pathogens in clinical samples.

[0032] (2) High distinguishability ensures the application of multiple internal standards: Through the second sequence-specific screening, it is ensured that the k-mers of each Spike-in sequence are not duplicated or the editing distance is large enough. After sequencing, the original sequence corresponding to each k-mer can be traced 100%, which solves the confusion problem when multiple internal standards are used simultaneously and adapts to the complex needs of clinical multi-gradient internal standard quality control, multi-target quantification, etc.

[0033] (3) High safety and stability: Arabidopsis thaliana is a non-pathogenic model organism with a stable and fully resolved genome sequence. It does not pose a risk of infecting humans, avoids the potential uncertainty of randomly synthesized sequences, and meets the safety requirements of clinical testing.

[0034] (4) Strong standardization and reproducibility: Key parameters such as sequence truncation length, k-mer size, Kraken2 conserved parameters, and edit distance threshold are clearly defined, forming a standardized process of "double screening + triple verification", which can repeatedly obtain qualified Spike-in sequences and is easy to prepare on a large scale (such as artificially synthesized plasmids and linear DNA fragments).

[0035] (5) Adaptable to all clinical testing scenarios: The physicochemical properties of the screened Spike-in sequence are adapted to the entire process of clinical metagenomic sequencing (extraction, amplification, sequencing), and can be used as a single internal standard or multi-gradient internal standard to achieve experimental quality control, absolute quantification and cross-batch data correction, thereby improving the accuracy and reliability of clinical pathogen detection. Attached Figure Description

[0036] Figure 1 This is a schematic diagram of the method flow of the present invention. Detailed Implementation

[0037] To make the technical solution of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0038] Example 1

[0039] Experimental materials and tools

[0040] Arabidopsis thaliana genome data: TAIR10 reference genome;

[0041] Sequence processing tools: seqkit v2.10.1, Python 3.12 (with Biopython library);

[0042] k-mer splitting tool: Python script;

[0043] Species classification tool: Kraken2 v2.1.6, paired with a custom microbial database (based on NCBI RefSeq 2025, containing 100,000+ genomes of common clinical pathogens).

[0044] Homology verification tools: BLASTn v2.14.0, NCBI NR / NT database (updated in 2024), Human Genome GRCh38;

[0045] Edit distance calculation tool: Python Levenshtein library (for calculating Hamming distance).

[0046] like Figure 1 As shown, the implementation steps of the clinical metagenomics Spike-in sequence design method are as follows:

[0047] Step 1: Candidate sequence library construction and download

[0048] Download the TAIR10 Arabidopsis thaliana whole genome sequence from the NCBI website (https: / / www.ncbi.nlm.nih.gov / datasets / genome / GCF_000001735.4 / ) and name it TAIR10.fasta; use the seqkit sliding command to segment each chromosome into 1000bp non-overlapping segments, as follows:

[0049] seqkit sliding -W 1000 -s 1000 -o candidate_sequences.fastaTAIR10.fasta

[0050] Approximately 119,665 1000bp candidate sequences were obtained, forming a candidate sequence library.

[0051] Step 2: k-mer splitting

[0052] The Python script is used to split each candidate sequence in candidate_sequences.fasta into 50bp k-mer segments, setting the k-mer length to 50, retaining all k-mer segments and associating them with their corresponding candidate sequence IDs. The command is as follows:

[0053] python get_kmers.py -i candidate_sequences.fasta -o 50mers.fasta -k50

[0054] By associating k-mers with candidate sequence IDs using Python scripts, the traceability of 951 k-mers corresponding to each candidate sequence is ensured.

[0055] Step 3: Screening for microbial homology with human sources (first screening)

[0056] (1) Kraken2 Database Construction: Download the genome sequences of clinical pathogenic microorganisms such as bacteria, fungi, viruses, and mycoplasma from NCBI RefSeq, and use the kraken2-build command to construct a custom microbial database. The command is as follows:

[0057] kraken2-build --db clinical_microbe_db --download-taxonomy

[0058] kraken2-build --db clinical_microbe_db --add-to-library refseq_microbe_fasta / *.fasta

[0059] kraken2-build --db clinical_microbe_db --build --threads 16.

[0060] (2) Kraken2 parameter settings: When running Kraken2, set confidence=0.4, as follows:

[0061] kraken2 --db clinical_microbe_db --confidence 0.4 --threads 16 --output 50mers.kraken --report 50mers.kreport 50mers.fasta.

[0062] (3) Sequence screening: The 50mers.kraken database was parsed using a Python script to screen candidate sequences for which none of the k-mers were classified as microorganisms or humans. If any k-mer of a candidate sequence was labeled by Kraken2 as "Bacteria", "Fungi", "Virus", "Human", or other microbial or human classifications, the sequence was directly filtered out. A total of 27,636 qualified candidate sequences were obtained after screening.

[0063] Step 4: Sequence-specific enhanced screening (secondary screening)

[0064] This embodiment uses "Strategy B (High Edit Distance Filtering)" to filter 10 sequences for multi-gradient internal standards:

[0065] (1) Extracting the k-mer set: Using a Python script, all 50bp k-mers of 119,665 sequences were extracted, forming a set containing approximately 119,665 × 951 ≈ 1.14 × 10 8 A global set of k-mers, recording the "k-mer - sequence ID" mapping.

[0066] (2) Greedy algorithm screening: ① Randomly select the sequence with sequence ID Seq_00123 as the reference sequence S1, and extract its 951 k-mers as the reference k-mer set K1; ② Traverse the remaining sequences, calculate the Hamming distance between the k-mer set of each sequence and all k-mers in K1, and select the sequence S2 with the minimum Hamming distance between all k-mers and the k-mers in K1 ≥ 3; ③ Using the k-mer set of S1 and S2 (K1∪K2) as the new reference, continue to screen the sequence S3 with the minimum Hamming distance between all k-mers in the set ≥ 10; ④ Repeat the above steps until 30 sequences S1-S3 are obtained, forming the second set of qualified sequences.

[0067] (3) Verification: Calculate the minimum k-mer distance between each pair of S1-S3. All of them are ≥10, which satisfies the distinguishability requirement.

[0068] Step 5: Spike-in sequence verification

[0069] (1) Secondary verification of microbial homology: S1-S3 were compared with microbial sequences in the NCBI NR / NT database using BLASTn. The E-value was 1e-5. All sequences showed homology of less than 65%, which met the requirements.

[0070] (2) Human homology verification: compared with the human genome GRCh38, the homology is less than 65%, and there is no risk of cross-reaction.

[0071] (3) Physicochemical property verification: The GC content of the 10 sequences is between 40% and 60%, which is suitable for PCR amplification and sequencing.

[0072] (4) Distinguishing verification: Sequencing simulation analysis showed that the k-mer traceability accuracy of each sequence reached 100%, with no confusion.

[0073] The embodiments described above are merely illustrative of several implementations of the present invention, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention. Therefore, the scope of protection of this patent should be determined by the appended claims.

Claims

1. A method for designing clinical metagenomic Spike-in sequences, characterized in that, Includes the following steps: S1. Construction of candidate sequence library: The Arabidopsis genome was selected as the sequence source. The whole Arabidopsis genome sequence was divided into candidate fragments of uniform length using a sequence truncation tool. The fragment length was set to 1000 bp to form a candidate sequence library. S2, k-mer splitting: For each 1000bp candidate sequence in the candidate sequence library, the k-mer splitting algorithm is used to split it into continuous short k-mer segments. The length of each k-mer is set to 50bp, forming a k-mer set corresponding to each candidate sequence. S3. Strict screening for microbial and human homology: Kraken2 software was used to annotate all k-mers of each candidate sequence for species classification. Conservative parameters were set to maximize the filtering of candidate sequences containing potential microbial or human homologous k-mers. S4. Sequence-specific enhanced screening: The final Spike-in sequence is selected from a set of qualified sequences that are strictly screened for microbial and human homology. S5, Spike-in sequence validation: Perform secondary validation on the qualified sequence set selected for enhanced sequence specificity to adapt it to the clinical metagenomic sequencing process.

2. The method for designing clinical metagenomic Spike-in sequences according to claim 1, characterized in that, In step S2, the k-mer splitting rule is: adjacent k-mers overlap by 49bp, so that each 1000bp candidate sequence can be split into 951 (1000-50+1) 50bp k-mers, forming the k-mer set corresponding to each candidate sequence.

3. The method for designing clinical metagenomic Spike-in sequences according to claim 1, characterized in that, In step S3, the specific screening process is as follows: (1) Kraken2 database configuration: Select the whole genome database and human genome database containing all common clinical pathogens such as bacteria, fungi, viruses, mycoplasma, chlamydia, and rickettsia from NCBI RefSeq, and update the database regularly to cover new pathogens; (2) Kraken2 conservative parameter settings: set the minimum k-mer matching length to 31bp, the minimum matching score to 0.95, turn off fast mode, and enable strict classification mode; (3) Screening criteria: For each candidate sequence, if any one of the 951 50bp k-mers is assigned to a microbial species or human source by Kraken2, the candidate sequence is identified as a potential interference sequence and filtered; only candidate sequences for which none of the k-mers are assigned to a microbial species or human source by Kraken2 are retained to form the first set of qualified sequences.

4. The method for designing clinical metagenomic Spike-in sequences according to claim 1, characterized in that, In step S4, the core requirement for screening is that there is no duplication or large edit distance between any two 50bp k-mers of each sequence, so that the original Spike-in sequence corresponding to each k-mer can be accurately traced after sequencing.

5. The method for designing clinical metagenomic Spike-in sequences according to claim 3, characterized in that, In step S4, the specific screening process is as follows: (1) Construct a global k-mer set: Extract the 50bp k-mer of all sequences in the first-stage qualified sequence set to form a global k-mer set, and record the sequence ID corresponding to each k-mer; (2) Selection of screening strategy: Strategy A: Traverse the global k-mer set and remove sequences containing duplicate k-mers; retain sequences whose k-mers are unique to the sequence and not shared by any other sequence to form a subset of k-mer sequences without duplicates; Strategy B: Perform pairwise comparisons on the qualified sequence set of the first screening, calculate the edit distance between all k-mers in the k-mer set corresponding to any two sequences, and screen the sequences that satisfy the condition that the minimum edit distance is ≥10 among all pairwise k-mer pairs of any two sequences; Choose one of strategy A and strategy B or use them in combination. (3) Determining the number of sequences: Based on clinical testing needs, select several sequences from the above subsets or screening results to form a second set of qualified sequences.

6. The method for designing clinical metagenomic Spike-in sequences according to claim 5, characterized in that, In step (2), if multiple sequences need to be selected, a greedy algorithm is used to select them step by step: first, randomly select one sequence as the benchmark, and then select sequences whose k-mer edit distances to the selected sequences are all ≥10, until the required number of sequences is reached.

7. The method for designing clinical metagenomic Spike-in sequences according to claim 1, characterized in that, In step S5, the specific process for sequence verification is as follows: (1) Secondary verification of microbial homology: The homology of the sequence with the microbial sequence in the NCBI NR / NT database was compared using the BLASTn tool, and the homology between the sequence and the microbial sequence was reduced to less than 70%; (2) Human sequence homology verification: BLASTn alignment with the human genome was performed to reduce homology to less than 50%; (3) Physicochemical property verification: Analyze the GC content of the sequence, retain the sequence with GC content between 40% and 60%, and adapt it to the PCR amplification and sequencing process; (4) Distinguishing verification: For the finally selected sequences, the k-mer repetition rate or edit distance between each pair is calculated again to verify that the distinguishability meets the standard.

Images