Primer design method of multiple PCR (Polymerase Chain Reaction) targeted sequencing technology

Abstract

The invention provides a primer design method of a multiplex PCR targeted sequencing technology, which designs a primer with high specificity and wide genome coverage rate by combining multiple indexes with fine screening of high-quality species genomes, effectively optimizes the formation of dimers of a primer system, and improves the specificity and efficiency of detection. According to the technology, the high specificity is guaranteed, meanwhile, the performance of a multiple primer system is remarkably improved, the common dimer problem is optimized, and the high efficiency of multiple pathogen detection is guaranteed.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, molecular biology, microbiology and sequencing technology, and in particular to a primer design method for multiplex PCR targeted sequencing technology. Background Art

[0002] With the development of science and technology, significant progress has been made in the fields of molecular biology, genetics, bioinformatics and sequencing technology. Among them, polymerase chain reaction (PCR), as one of the cornerstone technologies of molecular biology, amplifies specific DNA fragments in vitro, so that trace amounts of DNA samples can be detected; while the second and third generation sequencing technologies can sequence a large number of DNA fragments in parallel, thereby realizing rapid sequencing of the entire genome or transcriptome. Multiplex PCR (mPCR) technology is an extension of PCR technology, which allows multiple pairs of primers to be used simultaneously in a PCR system for simultaneous amplification of multiple target sequences. Compared with traditional PCR, this method has the advantages of high efficiency, high throughput and low cost, and can meet the clinical needs for simultaneous detection of multiple pathogens and microorganisms, thereby greatly improving diagnostic efficiency and accuracy. However, in practical applications, these technologies still have some challenges and limitations, which are mainly reflected in the following points: 1. The problem of genome library quality screening: The construction of a high-quality species genome library is the foundation for the subsequent selection of species-specific segments and the design of ideal primers. Existing genome screening methods mostly rely on genome assembly-level tags (such as representative genome, reference genome, complete genome, chromosome, etc.), lacking more detailed screening criteria (such as genome sequencing time, host, origin, completeness, etc.), resulting in uneven genome quality in the genome library.

[0003] 2. Primer compatibility issues: When designing multiplex PCR primers to detect target pathogens, since the genomes of target pathogens in clinical samples are diverse, the primers themselves need to have high compatibility, that is, the primers for each target pathogen can copy the diverse genomes with slightly different sequences under the pathogen. Existing methods usually first look for highly conserved species-specific fragments in the genome of the species, and then use the specific fragments as templates to design primers. However, for species with more diverse genomes, the process of finding specific fragments is often not smooth, and it is necessary to continuously adjust the genome library for multiple rounds of specific sequence searches. The workload is large and it is difficult to obtain specific fragments with a fragment length that can meet the primer design requirements and high coverage.

[0004] 3. Primer dimer problem: In multiplex PCR system, due to the large number of primers, the dimer binding between primers will reduce the efficiency of the system. Common methods only screen and optimize primer dimers based on the binding free energy value, without fully considering the influence of factors such as the base length of the primer dimer binding region, the distance between the starting and ending points of the dimer binding region and the 3' end of the two primers on the dimer binding. Summary of the invention

[0005] Purpose of the invention: The purpose of the present invention is to provide a primer design method for multiplex PCR targeted sequencing of pathogens, so as to avoid a lot of duplication of work in obtaining species-specific sequences with diverse intraspecies genomes and provide ideal compatibility.

[0006] Technical solution: The present invention provides a primer design method for multiplex PCR targeted sequencing technology, which comprises the following steps in sequence: S100: Determine the best reference genome and construct a high-quality genome library of the target pathogen; S200: constructing a phylogenetic tree for the high-quality genome library; S300: traversing the optimal reference genome with a fixed step length to obtain candidate primer pairs; S400: Evaluate the candidate primer pairs, including scoring of coverage, specificity and dimer binding.

[0007] The target pathogens described in the present invention refer to a variety of organisms targeted by pre-designed primers in a targeted multiplex PCR primer system, including but not limited to viruses, bacteria and fungi.

[0008] The primer design and optimization method provided by the present invention is applicable to the field of high-throughput molecular diagnosis, and is particularly suitable for multiplex PCR to detect multiple target pathogen targets through a single PCR reaction. Due to the diversity of target pathogen genomes in clinical samples, it is particularly important to select primers with high compatibility to ensure that diverse genomes with slight sequence differences can be effectively copied. Building a high-quality and diverse species genome library is the basis for ensuring high compatibility of primers. In this way, while ensuring that the primers highly cover the genome library, it can be inferred that the primers also have high compatibility with various genomes in the real world. At the same time, when designing primers, it is also necessary to improve the compatibility of primers to diverse genomes and adopt efficient design strategies to reduce cumbersome workloads. The present invention can solve these two difficulties well.

[0009] The present invention constructs a high-quality genome library of target pathogens, which is different from the prior art that is screened according to the order of genome assembly levels, taking into account multiple criteria such as genome sequencing time, host, place of origin, and integrity, so that when screening, the genome with general sequencing quality produced due to early sequencing year and insufficient sequencing technology can be filtered out, or the genome of the host and place of origin and the Chinese population is too far from the source. If these genomes are present in the genome library, the evaluation of primer compatibility will be interfered. At the same time, in order to ensure that the genome library has diversity, the present invention also uses the method of constructing a phylogenetic tree to cluster the genome library on kinship, ensuring that the genome in the library has neither most of the genome concentrated in a certain subspecies branch, poor diversity, nor individual distant evolutionary branches missing when designing primers. The construction of a high-quality genome library is based on the best reference genome, and other genomes of the same species are obtained according to assembly level, genome sequencing time, host, place of origin, integrity, etc., and the sequencing time of the same species genome is required to be close to the current year, and there is a certain number of samples.

[0010] In the primer design stage, in order to improve the compatibility of primers to the diverse genomes of target species, the common method in the industry is to first obtain the conserved fragments shared by multiple genomes in the genome library, and then design primers on these conserved fragments. In order to obtain these highly conserved fragments, multiple rounds of genome sequence alignment and clustering are required to find the common fragments. This process is cumbersome and time-consuming, especially for some species with high genome diversity, after a lot of work, it is still not possible to successfully find sequences with both high conservatism and high species specificity. The present invention first determines an optimal reference genome in the species genome library, and then clusters the diverse species genome library with a phylogenetic tree. If there is a cluster of distant sources, it is grouped separately, and an optimal reference genome is also confirmed in the grouping. Design candidate primers on the optimal reference genome, and then calculate the coverage of these candidate primers in the genome library. For the distant genome grouping, design candidate primers on the reference genome on the grouping, and then calculate the genome coverage of the candidate primers in the grouping. This method of first grouping, then designing primers, and finally calculating primer coverage reduces the large amount of long fragment alignment work when obtaining common fragments. It is more concise and also ensures high primer coverage.

[0011] Furthermore, the step S100 specifically includes the following steps: S110: Obtain the best reference genome of the target pathogen; S120: obtaining a number of genomes of other species with the most recent sequencing time, wherein the genomes of other species are obtained in order of priority from high to low assembly levels; S130: Screening out other genomes of the same species that have high completeness and meet the sample quantity requirements. When the sample quantity is insufficient, preferentially extending the sequencing time of step S120.

[0012] Among them, the optimal reference genome is determined based on multiple dimensions such as genome source, genome assembly level, sequencing time (year) and CheckM value to ensure that the screened genome is of high quality and representative of the species.

[0013] Specifically, the method for determining the best reference genome is: selecting a best reference genome from the RefSeq database or the GenBank database in the priority order of Representative Genome>Reference Genome>Complete Genome>Chromosome>Scaffold>Contig; When there are multiple genomes or genome structures at the same level, the genome or genome structure with the most recent sequencing year and the highest CheckM value is selected; Among them, the priority of sequencing year is greater than CheckM value; and if the sequencing years and CheckM values ​​of multiple genomes or genome structures of the same level are the same, one is selected arbitrarily.

[0014] Specifically, species in the RefSeq database usually contain at least one genome of Representative Genome and Reference Genome. If the number of Representative Genome or Reference Genome of a species is greater than 1, the genome with the closest sequencing year and the highest CheckM value is selected as the best reference genome, and the sequencing year has a higher priority than the CheckM value; if the sequencing year and CheckM value are the same, one is selected at random. If the species does not have a Representative Genome or Reference Genome in the Refseq database, the best reference genome is confirmed in the GenBank database in the same way as above. If there is no Representative Genome and Reference Genome in both the Refseq and GenBank databases, the genome structure is screened downward in the order of Complete Genome>Chromosome>Scaffold>Contig. When the number of genomes / genome structures at the same level is greater than 1, the genome with the closest sequencing year and the highest CheckM value is selected as the best reference genome, and the sequencing year has a higher priority than the CheckM value; if the sequencing year and CheckM value are the same, one is selected at random.

[0015] The CheckM value used to confirm the optimal reference genome in the present invention is an indicator for measuring the completeness of the genome.

[0016] Among them, for pathogens that are taxonomically classified as "species", the number of samples shall not be less than 10; for pathogens that are taxonomically classified as "species complex" or "genus", the number of samples shall be 50.

[0017] It should be noted that step S120 gives priority to obtaining genome samples with complete assembly level. When the number of samples is insufficient, the sequencing time of step S120 is preferentially extended to obtain more samples; when the extended sequencing time still cannot obtain sufficient samples, chromosome, scaffold, contig and other genes are downloaded according to the assembly level priority until the number of samples reaches the standard.

[0018] As a further optimization of the present invention, each target pathogen in the high-quality genome library is tested for integrity using CheckM2 software to screen genomes with integrity greater than or equal to 95%.

[0019] In particular, the genomes of some species are more complex and CheckM2 cannot be scored, and the CheckM items of NCBIRefSeq / GenBank are all 0. In this case, this indicator can be ignored.

[0020] The phylogenetic tree described in the present invention is also called a molecular evolutionary tree. It is a visualized tree diagram constructed by comparing the numerical values ​​of the sequence differences of biological macromolecules (such as proteins and nucleic acids). It describes the kinship between different objects at the sequence level and is a bioinformatics analysis method. In the prior art similar to the subject matter of the present application, there is little disclosure of constructing a phylogenetic tree to test the evolutionary kinship of each genome in the high-quality genome of a species, and applying it in the primer design method of multiple PCR targeted sequencing technology. The present invention is based on the phylogenetic tree to "cluster" the genome from the perspective of genome sequence differences to reflect the kinship between genomes. Specifically, the present invention uses Parsnp software to construct a phylogenetic tree, which first aligns the maximum similar sequence parts of each genome in the species to form a "core alignment region", and then detects single nucleotide polymorphisms in the "core alignment region". These single nucleotide polymorphisms are usually regarded as the most reliable mutations for studying evolutionary relationships; then the maximum likelihood algorithm is used to construct an evolutionary tree based on the differences in the "core alignment region" to reflect the kinship of these genomes.

[0021] Further, the genome evolution tree obtained in step S200 meets the following requirements: the branches of the genome evolution tree are uniform, there are no isolated branches with obviously long evolutionary distances, and there is no situation where most genomes are densely distributed in the same branch. In the phylogenetic tree, the leaf node represents the terminal node of the phylogenetic tree, and there is no further branch at this node. In this study, each leaf node represents a genome participating in the construction of the evolutionary tree. The initial evolutionary tree first analyzes the similarities and differences of multiple genomes with similar sequences to form an independent small branch structure, and multiple small branch structures then form branch extensions tending to the root node. The branch length represents the relationship between two nodes. The longer the branch, the farther the relationship. The uniform branching of the evolutionary tree is determined as follows: in the phylogenetic tree, in the small branch structure with the most leaf nodes, the number of leaf nodes is less than 30%~50% of the total number of genomes. Furthermore, if the ratio of the length of some evolutionary branches to the average length of all evolutionary branches is greater than 10-100 times, the genomes on these branches are considered distant genomes. If the number of genomes on the branches is > 0.01-10% of the total number of genomes, and the quality of these genomes is good, these genomes are grouped according to the evolutionary tree branches, and primers are subsequently designed within the group. If the number of genomes on the branches is > 0.01-10% of the total number of genomes, these genomes are directly deleted.

[0022] As a preferred embodiment of the present invention, whether the genomes of the obviously distant isolated branches are removed from the high-quality genome library, and whether or not the genomes with obviously distant evolutionary distances exceeding a certain proportion of the total number of genomes are grouped according to the evolutionary tree branches is evaluated based on 5% of the total number of genomes. It should be understood that in the design of multiplex PCR primers, in order to cope with the diversity of the target pathogen genomes in clinical samples, it is necessary to ensure that the primers have high compatibility and can cover the diverse genome sequences of each target pathogen. Therefore, it is necessary to select genomes with high sequencing quality and universality to construct a high-quality genome library. This requires that the selected genomes should not be concentrated in a single evolutionary lineage (such as a certain subspecies) to avoid insufficient universality. Similarly, when a group of genomes is too distant in evolutionary relationship and the number exceeds 5% of the total number, it should be divided into several groups according to the evolutionary relationship, and primers should be designed to cover them separately.

[0023] For step S300, in step S300, 300 bp is used as the sliding window length and 20 bp is used as the sliding step length, and the whole genome is sliced ​​in a shingled manner, and primers are designed on each fragment using primer3 software.

[0024] Preferably, for the second generation sequencing technology, the amplicon length is 180-300 bp, T m The value range is 55-65°C, the primer arm length is 18-25bp, and the primer arm GC base content is 40-60%.

[0025] For the third-generation sequencing technology, those skilled in the art may make adaptive changes to the design of the sliding window length and the sliding step size based on the above-mentioned second-generation sequencing technology without departing from the inventive principles of the present application, and all of these are within the protection scope of the present application.

[0026] For step S400, the coverage evaluation method is to match the candidate primer pairs to each other genome of the same species in the high-quality genome library, and calculate the primer coverage: Primer coverage = number of qualified matching genomes / total number of genomes × 100%.

[0027] Among them, a qualified match at least meets the following conditions simultaneously: both the upstream primer and the downstream primer can achieve a match, the matching length of the primer amplicon is not less than 95% of the total length of the amplicon, and the primer amplicon length is 100 - 300 bp.

[0028] The specific evaluation method comprises the following steps: S421: aligning the candidate primer pairs in pairs in the NCBI nt library against all genomes in the library; S422: After alignment, a primer insert fragment is obtained, wherein the primer insert fragment refers to a nucleic acid sequence defined and amplified by a pair of complementary oligonucleotide primers during PCR, wherein the starting end of the primer insert fragment is defined by the 3' end of the forward primer, and the ending end is defined by the 5' end of the reverse primer; S423: inputting the primer insert fragment into Kraken 2 software for verification, thereby screening primer pairs with qualified specificity; S424: When the species of the genome to which the primer insert fragment belongs is consistent with the primer target species, the species discrimination result of Kraken 2 must be that species; or when the species of the genome to which the primer insert fragment belongs is inconsistent with the primer target species, the species discrimination result of Kraken 2 must be consistent with the genome species.

[0029] The above method allows primers to bind to sequences of non-target species, while ensuring that non-target species can be correctly identified based on primer inserts in Kraken 2. This relaxes the requirements for primer specificity and does not generate false positives due to species misreporting, thus improving primer utilization.

[0030] In a multiple primer system, one of the common reasons that affect the amplification efficiency of the system is the dimer binding between primers, that is, the primers in the system bind to each other due to base complementarity, which makes the primers unable to bind to the sequence of the pathogen, reducing the amplification efficiency of the primers. The common method in the industry is to first enumerate all the dimers that may produce binding in the system and calculate the free energy of these dimer binding. When the free energy is outside the safety threshold range (<- 4 kcal / mol), it is considered that the two primers are more likely to combine to produce dimers in real experiments. Or use some global optimal algorithms, such as greedy algorithms, to combine primers, try to select new primers that have the least binding with primer dimers in the existing system, and cycle multiple times until the number of primers reaches the standard. The calculation of dimer binding free energy involves base pairing, stacking interaction, ion concentration, and temperature indicators. It only measures the ease of binding of this segment of bases at a specific temperature and solution concentration, but ignores that the binding of the dimer is also related to the distance between the dimer region and the 3' end of the two primers. The closer the dimer binding region is to the 3' end, the worse the effect. The global optimal algorithm often results in a fixed set of primer combinations. When some primers in the system need to be replaced, the situation will be very troublesome, and the algorithm may need to be re-run.

[0031] The present invention proposes a new primer dimer evaluation method, which not only takes into account the dimer binding free energy that is of general concern to those skilled in the art, but also integrates factors such as the length of the primer dimer binding region, the number of GC bases in the binding region, and the distance between the primer dimer binding site and the 3' end of each primer. Traditional methods often rely on algorithms to evaluate primer pairs one by one in sequence, resulting in a lack of flexibility and controllability in the results. The evaluation method provided by the present invention is more intuitive, allowing researchers to freely select primer pairs, rather than simply relying on a fixed combination automatically generated by an algorithm, thereby providing greater flexibility and freedom of choice.

[0032] Specifically, the dimer binding score is: ; When any one of the candidate primer pairs No less than 10 3 When , the primer pair corresponding to the primer is screened out; "Dimer binding region" refers to the double-stranded region formed by partial or complete complementarity between two complementary oligonucleotide primers during the polymerase chain reaction (PCR) process. This region starts from the non-complementary region at the 3' end of any primer and may extend to the 3' end of the other primer or any other position, terminating at the end point of complementary pairing. In this segment, individual base mismatches are allowed, but sufficient complementarity is maintained as a whole to form a stable double-stranded structure. GC-count represents the number of GC bases in the dimer binding region (mismatched bases are not counted), which is based on the first complementary binding base from the 5' end to the 3' end of one of the primers as the starting point, and the last complementary binding base as the end point. The start-end region can accept up to 2 base mismatches; L represents the base length from the start to the end of the dimer binding region; d1 represents the distance from the start of the dimer binding region to the 3' end of the nearest primer, and d2 represents the distance from the end of the dimer binding region to the 3' end of the nearest primer.

[0033] The above indicators are jointly involved in dimer evaluation, which can more comprehensively evaluate the difficulty of dimer binding area and the impact on the entire system. In the end, each candidate primer will have a score, which is the sum of the scores of all possible dimers generated by this primer and other primers in the system. The higher the score, the more serious the dimer binding situation. Researchers can directly screen and combine according to the score, which is more flexible than conventional evaluation algorithms.

[0034] Based on the dimer evaluation method provided by the present invention, when the total score of any primer in a pair of primers is greater than or equal to 10 3 , then the pair of primers should be considered to be removed from the system. For example, the F-Primer and R-Primer shown below both have a dimer score of 4096 using the dimer provided by the present invention, and should be evaluated as primer pairs that need to be eliminated.

[0035] Score: 6, Delta G = - 8.17 kcal / mol F-Primer: CGTTCAGTACAATG CGGCCG R-Primer: GCCGGC GTAACATGACTTGC The underlined part of the above primer is the dimer junction region, the binding length is 6 bases, and all of them are GC bases. The dimer binding free energy is -8.17 kcal / mol, which is less than the threshold value of -4 kcal / mol generally specified by those skilled in the art, and belongs to a dimer that is easy to bind in the real world. It can be seen that the dimer evaluation method provided by the present invention can objectively evaluate the dimer binding situation and use it to screen primers.

[0036] It should be noted that, in the technical solution claimed in the present application, the order of primer coverage, specificity and dimer binding score is used as the preferred evaluation method, that is, the primer pairs screened by the first evaluation are used for the subsequent evaluation. Nevertheless, those skilled in the art can change the evaluation order of the above content at will if necessary, which is within the scope of protection of the present application.

[0037] The primer design method provided by the present invention includes determining the best reference genome of the target pathogen, building a high-quality genome library of the pathogen, designing primers, and evaluating the coverage, specificity and dimer binding of the primer pairs. The method first traverses the best reference genome of the target pathogen to design primer pairs, and then aligns the primers to the gene library to calculate the coverage, thereby avoiding a lot of duplication of work when obtaining species-specific sequences with more diverse genomes within the species; then the specificity of the candidate primer pairs is tested to ensure that the fragments copied by the primer pairs in the real sample can correctly identify the species to which they belong; finally, the candidate primer pairs are scored using a self-developed dimer scoring algorithm to reduce dimer binding in the primer system and improve the working efficiency of the primer system. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 Flowchart for designing primer system for multiplex PCR targeted pathogen detection; Figure 2 This is the evolutionary tree of the high-quality genome of Staphylococcus worteri; Figure 3 The evolutionary tree of the high-quality genome of Candida tropicalis; Figure 4 Phylogenetic tree of the high-quality measles virus genome. DETAILED DESCRIPTION

[0039] In order to make the technical solution of the present invention clearer, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0040] Example 1

[0041] This example uses Staphylococcus worderi (bacteria), Candida tropicalis (fungus) and measles virus (virus) as target pathogens and is described in detail based on Example 1.

[0042] like Figure 1 As shown, the present invention provides a primer design method for multiplex PCR targeted sequencing technology, and the method comprises the following steps in sequence: S100: Determine the best reference genome and construct a high-quality genome library of the target pathogen; S200: constructing a phylogenetic tree for the high-quality genome library; S300: traversing the optimal reference genome with a fixed step length to obtain candidate primer pairs; S400: Evaluate the candidate primer pairs, including scoring of coverage, specificity and dimer binding.

[0043] S100. Construction of high-quality genomic libraries

[0044] S110. Obtain the best reference genome for the target pathogen

[0045] Obtain the assembly_summary_refseq.txt file from the RefSeq database of NCBI. Confirm the target species based on the species name field (organism_name) combined with the species ID field (species_taxid), and confirm the best reference genome based on the genome category field (refseq_category) as "reference" or "representative", and complete the download of the best reference genome for each target species.

[0046] The species name of Staphylococcus worderi is Staphylococcus warneri , species ID: 1292; the species name of Candida tropicalis is Candida tropicalis , species ID: 5482; the species name of measles virus is Measles morbillivirus , species ID: 11234. All three species can be directly confirmed and downloaded to representative genome from the RefSeq database.

[0047] S120. Obtain several other genomes of the same species that were sequenced recently

[0048] The target pathogen is determined according to the species name field (organism_name) obtained in step S110 in combination with the species ID field (species_taxid); and through the NCBI information page of the genome, other genomes of the same species whose sequencing information was published later than 2018 and whose host field (host) is "Homo sapiens" or "Human" are screened.

[0049] Among them, as a preferred embodiment, source fields (source) including "soil", "water" and other genomes that are not related to humans or mammals can be screened out. These genomic information is of little use in genomic research on human infectious pathogens.

[0050] Download other genomes of the same species in the order of assembly level complete genome>chromosome>scaffold>contig from high to low priority. When the target pathogen is at the "species" level, the total number of genomes must be no less than 10; when the target pathogen is at the "species complex" or "genus" level, the total number of genomes must be >= 50.

[0051] If the number of genomes at the complete genome level of the current pathogen still does not meet the minimum number standard (10) after downloading all, the sequencing time standard will be relaxed by 2 years (for example, "2018 to date" will be relaxed to "2016 to date" for the first time). If the number of genomes has not yet reached the standard, they will be downloaded in the order of assembly level priority until the number of genomes reaches the standard.

[0052] S130. Screen out genomes of other species with high completeness and sufficient sample size

[0053] The best reference genomes of the three target pathogens were tested for genome integrity using CheckM2 software (CheckM2 v1.0.2). The test standard was that the score of the integrity index CheckM was not less than 95%.

[0054] If the number of unqualified genomes of the current pathogen is not less than 15%, the candidate genomes are added according to S120 and the completeness is performed. In this embodiment, because the species level of the three target pathogens is species, the total number of genomes of each species is required to be no less than 10. The genome of Staphylococcus worderi can be scored using CheckM2; the genomes of Candida tropicalis and measles virus cannot be scored using CheckM2, and the CheckM in the NCBI RefSeq and GenBank databases are both 0, so this indicator is not considered for these two pathogens. Table 1, Table 2, and Table 3 are the genome information and CheckM scores of Staphylococcus worderi, Candida tropicalis, and measles virus, respectively.

[0055] Table 1 List of high-quality genome libraries of Staphylococcus worteri

[0056] Table 2 List of high-quality genome libraries of Candida tropicalis

[0057] Table 3 List of high-quality genome libraries of measles virus

[0058] S200. Use Parsnp to construct a genome evolution tree for high-quality genome libraries For each genome screened by the target pathogen, an evolutionary tree was constructed using Parsnp software (Parsnp v1.7.4), and then visualized using the tools provided on the online website to test whether the genome was evenly distributed from the perspective of evolutionary relationships.

[0059] The best reference genome of each target pathogen is used as a reference to construct a genome evolutionary tree with all other genomes of the same species. The branches of the evolutionary tree should be relatively uniform, with no isolated branches with significantly long evolutionary distances, and no dense concentration of all genomes on the same branch. Specific judgment method: In the small branch structure with the most leaf nodes, the number of leaf nodes should be less than 50% of the total number of genomes. For branches with a ratio of branch length to average branch length greater than 100, if the leaf nodes contained are greater than 10% of the total number of genomes, and the quality of the genomes is good, these genomes are grouped according to the evolutionary tree branches, and primer design within the group is subsequently performed; if the number of leaf nodes contained is less than or equal to 10% of the total number of genomes, these genomes are directly deleted.

[0060] The phylogenetic trees of Staphylococcus wortii, Candida tropicalis, and measles virus are shown in Figure 2 , Figure 3 , Figure 4 . Figure 2 In the figure, the dotted box circles a small branch structure with the most leaf nodes, with a total of 5 leaf nodes, which does not exceed 50% of the total number of genomes. The average length of the evolutionary tree branches is 0.114, and the longest branch length is 0.389, which does not exceed 100 times the average length. The results show that the genome evolutionary tree is evenly distributed and there are no obviously isolated genomes, which is considered to be qualified for the construction of the genome library. Similarly, it can also be judged that Figure 3 , Figure 4 The evolutionary trees are all in line with expectations. The high-quality genome libraries of the three species are qualified, and there is no need to delete genomes or group certain genomes to design primers. It should be noted that some branches of the phylogenetic tree are not marked with branch lengths because the software calculates them. They are not included in the calculation of the longest branch length and the average branch length, and are not deleted from the gene library. They are included in the calculation when calculating the coverage of the primers on the genome in step S400, thereby measuring the coverage of the primers on these genomes.

[0061] S300. Traverse the optimal reference genome with a fixed step size to obtain candidate primer pairs Candidate primer pairs were designed iteratively on the best reference genome for each target pathogen.

[0062] Select the best reference genome in each target pathogen library. This example takes the second-generation sequencing application scenario as an example: with 300 bp as the sliding window length and 20 bp as the sliding step length, perform shingled slicing on the whole genome, and use primer3 software to design primers on each fragment. If targeted amplification is performed on the second-generation sequencing platform, the primer design parameters are: the full length range of the amplicon is 180~300bp; T m The value range is 55-65℃, where T m The value was predicted by the calculate_tm() function of the software primer3; the primer arm length range was 18-25bp; the primer arm GC base content range was 40-60%, and the primer arm did not contain >= 5 consecutive identical bases. After design, the primer pairs that were exactly the same were removed.

[0063] When replaced with a third-generation sequencing platform, the present invention can match the requirements by changing the sliding window, sliding step length and primer design parameters.

[0064] S400. Evaluating the candidate primer pairs, including coverage, specificity, and dimer binding scores S410. Coverage evaluation Calculate the coverage of the candidate primer pair of each species in the high-quality genome library of the species. Match the candidate primer pair obtained in step S300 to each other genome of the same species in the high-quality genome library, and calculate the primer coverage: Primer coverage = number of qualified matching genomes / total number of genomes × 100%.

[0065] The criteria for qualified matching are: both the upstream primer and the downstream primer can match the genome; the primer amplicon length is between 100-300bp; and the matching length of the amplicon is ≥ the total length of the amplicon × 95%; at least 2 consecutive bases at the 3' end of the primer arm completely match the genome.

[0066] After coverage evaluation, the primer information of the three target species that met various criteria is shown in Tables 3, 4 and 5. These candidate primers were involved in the subsequent optimization steps.

[0067] Table 4 Primer information with qualified coverage of Staphylococcus worderi

[0068] Table 5 Primer information with qualified coverage of Candida tropicalis

[0069] Table 6 Primer information with qualified measles virus coverage

[0070] S420. Specificity evaluation The specificity evaluation includes the following steps: S421: comparing the candidate primer pairs in pairs in the NCBI nt library to all genomes in the library, i.e., primer-blast; S422: After alignment, a primer insert fragment is obtained, wherein the primer insert fragment refers to a nucleic acid sequence defined and amplified by a pair of complementary oligonucleotide primers during PCR, wherein the starting end of the primer insert fragment is defined by the 3' end of the forward primer, and the ending end is defined by the 5' end of the reverse primer; S423: inputting the primer insert fragment into Kraken 2 software for verification, thereby screening primer pairs with qualified specificity; S424: When the species of the genome to which the primer insert fragment belongs is consistent with the primer target species, the species discrimination result of Kraken 2 must be that species; or when the species of the genome to which the primer insert fragment belongs is inconsistent with the primer target species, the species discrimination result of Kraken 2 must be consistent with the genome species.

[0071] Species identification was performed using the software Kraken 2 (Kraken2 v2.1.3). It was necessary to add a custom high-quality genome of the species to the library of Kraken 2 in advance and mix it with the standard species library that comes with Kraken 2 to build a library. The first 42bp of the candidate primer insert fragment was used for species identification in Kraken 2. When the species of the genome to which the insert fragment belongs is consistent with the primer target species or is a child of the target species, the species discrimination result of Kraken 2 must also be the target species or a child of the target species; or if the species of the genome to which the insert fragment belongs is inconsistent with the primer target species and is not a child of the target species, the species discrimination result of Kraken 2 must be consistent with the genome species. In this way, the primer pair is allowed to have copies of non-target species, but the species must be correctly identified in Kraken 2 without generating false positives.

[0072] Tables 7 to 9 show the proportions of the four situations of “copying the target species and reporting correctly”, “copying the non-target species and reporting correctly”, “copying the target species and reporting incorrectly”, and “copying the non-target species and reporting incorrectly” after the candidate primers were simulated amplified in the NCBI nt library and then identified the species. The candidates were considered qualified when the proportion of “copying the target species and reporting correctly” was greater than or equal to 90%, and the proportions of “copying the target species and reporting incorrectly” and “copying the non-target species and reporting incorrectly” were both 0, or there were special circumstances that could be explained.

[0073] Among the candidate primers for Staphylococcus vortex in Table 7, SW_1_F / SW_1_R and SW_3_F / SW_3_R were considered to be unqualified primers with low efficiency because they "copied non-target species and reported incorrectly" and the proportion of "copied target species and reported correctly" was less than 90%, and the remaining primers were qualified; all candidate primers for Candida tropicalis in Table 8 were qualified in specificity; it should be particularly noted that among the candidate primers for measles virus in Table 9, although each pair of primers had a small proportion of "copied non-target species and reported incorrectly", the reason was that the artificial synthetic plasmid with taxid 32630 was amplified and reported as the target species, or the artificial synthetic vector with taxid 2877938 was amplified and reported as the target species. Such artificial synthetic species do not exist in real-world human samples, and therefore do not affect the amplification and species reporting of the primers, so the candidate primers are considered to be qualified.

[0074] Table 7 Specificity verification information of candidate primers for Staphylococcus worteri

[0075] Table 8 Specificity verification information of candidate primers for Candida tropicalis

[0076] Table 9 Specificity verification information of candidate measles virus primers

[0077] S430. Dimer binding score All qualified candidate primers in step S420 are enumerated using the software MFEprimer (MFEprimer v3.2.7) to enumerate all possible dimer binding situations. The dimer binding length must be greater than or equal to 3 bases and the binding region is allowed to have a mismatch of less than or equal to 2 bases.

[0078] On this basis, each dimer binding is scored, and the scoring rules are: ; Wherein, GC-count indicates the number of GC bases in the dimer binding region, L indicates the base length of the dimer binding region; d1 indicates the distance from the starting point of the dimer binding region to the 3' end of the nearest primer, and d2 indicates the distance from the end point of the dimer binding region to the 3' end of the nearest primer.

[0079] Scoring of a primer pair Equal to two primers All primer pairs are based on The values ​​are sorted in ascending order, The smaller the value, the less dimer binding of the primer pair. In the multiplex PCR system, the present invention considers that the ≥10 3 is unqualified.

[0080] Tables 10, 11 and 12 are the dimer scores of all qualified candidate primers in step 7, among which CT_1_R and MM_2_F have dimer binding and the score is 8192 greater than 10 3 , which is considered to be a dimer that is easier to form, so the two pairs of primers CT_1 and MM_2 are eliminated in pairs, and the remaining candidate primers are qualified. The dimer formed by CT_1 and MM_2 is shown below, where the two bases marked in brackets represent the mismatched bases in the dimer binding region.

[0081] CT_1_R:GTAGAATAGTTAACT CGG(GA)TCTGT MM_2_F: GCC(AG)AGACA GTAATTGGAGGGG Table 10 Dimer verification information of candidate primers for Staphylococcus wortelii

[0082] Table 11 Dimer verification information of candidate primers for Candida tropicalis

[0083] Table 12 Dimer verification information of candidate primers for measles virus

[0084] Example 2 In this example, the primers screened in Example 1 were subjected to a single-pair PCR wet test to verify the feasibility of the method described in Example 1.

[0085] (1) Experimental reagents and materials Human cerebrospinal fluid samples; sample pathogens; self-developed nucleic acid extraction kit; reverse transcription kit (Abclonal, catalog number RM20428); SYBR Green QPCR kit (Aikuri Biotechnology, catalog number AG11701).

[0086] (2) Test equipment TL2010S medium-throughput tissue grinder (Dinghaoyuan Technology, catalog number 0401265); automatic nucleic acid extractor (hema, catalog number E96-II); LightCycler 480Ⅱ real-time fluorescence quantitative PCR instrument (Roche, catalog number LC480); PCR instrument (Longji, catalog number A300); vortex mixer vortex 5 (Qilin Bell, catalog number HD201805135); MiniL-12G large-capacity mini centrifuge (Maigao, catalog number MinilL-12G), etc.

[0087] (3) Primer synthesis and preparation In the single-plex PCR test, a universal sequence needs to be added to the specific primer sequence during primer synthesis, wherein the upstream primer is composed of the universal sequence, UMI and the specific upstream primer, and the downstream primer is composed of the universal sequence and the specific downstream primer. The primers are synthesized by Sangon Biotechnology Co., Ltd. and Acry Biotechnology Co., Ltd.

[0088] (4) PCR reaction system and reaction procedure The amplification efficiency of the primers was evaluated using the QPCR system. Using pure bacteria or plasmids as templates, two groups, a positive template group and a blank control group, were set up. The QPCR reaction was performed according to the QPCR reaction system in Table 13 and the QPCR reaction procedure in Table 14. The amplification efficiency of a single pair of primers was evaluated based on the cycle threshold (Ct) and melting temperature (Tm value). In particular, because the measles virus is an RNA virus, the sample used is RNA, which needs to be reverse transcribed before the QPCR reaction. X in the table indicates that the specific input amount is determined according to the actual situation.

[0089] Table 13 PCR reaction system

[0090] Table 14 PCR reaction program

[0091] (5) Qualification criteria for single-plex PCR primers The amplification efficiency of a single pair of primers was evaluated using SYBR Green dye fluorescence quantitative PCR, which was divided into excellent, good, and poor. The evaluation factors were: ① Blank template Ct-positive template Ct≥5; ② No amplification of blank template or blank template Ct≥25; ③ | Positive template T m Value - Blank Template T m Value|≥2.

[0092] The results of single-plex PCR with all primers are shown in Table 15. ...

[0093] Table 15 Single-plex QPCR test results of target pathogen primers

[0094] Table 15 (continued) Single-plex QPCR test results of target pathogen primers

[0095] A large number of experiments have shown that primers (pairs) with excellent or good amplification efficiency can play a role in the multiplex PCR system and can effectively amplify the target pathogen. Therefore, the single pair PCR verification of the primers screened in Example 1 is also qualified.

[0096] Example 3 Wet-lab validation of the multiple primer system in real clinical samples

[0097] To further verify the effectiveness of the design and optimization of the multiple primer system of the present invention, the qualified primers of the three target pathogens mentioned above were further subjected to wet test verification of the mixed system in real clinical samples. Due to the complexity of real clinical cerebrospinal fluid samples and the generally low concentration of pathogens, the wet test verification of this example can more fully verify the high efficiency and high sensitivity of the primers designed by the present invention.

[0098] (1) Experimental reagents and materials Positive human cerebrospinal fluid samples known to contain the target pathogen; self-developed nucleic acid extraction kit; reverse transcription kit (Abclonal, Catalog No. RM20428); first-round amplification kit (TaKaRa, Catalog No. RR062B); second-round amplification kit (Vazyme, Catalog No. N616-02); dsDNA quantification kit (Yeason, Catalog No. 12640ES76); ssDNA quantification kit (Yeason, Catalog No. 12645ES76); Hieff NGS DNA Selection Beads (Yeason, 12601ES75); DNBSEQ-E25RS High-throughput Sequencing Reagent Set (FCL SE100) (MGI, Catalog No. 940-000573-00); UltraPureDNase / RNase-Free Distilled Water (Invitrogen, Catalog No. 10977015), etc.

[0099] (2) Test equipment TL2010S medium-throughput tissue grinder (Dinghaoyuan Technology, catalog number 0401265); automatic nucleic acid extractor (hema, catalog number E96-II); PCR instrument (Longji, catalog number A300); Qubit 3.0 fluorescence quantitative instrument (Thermo Fisher, catalog number Q33216); DNBSEQ-E25RS gene sequencer (MGI, catalog number 900-000473-00); vortex mixer vortex 5 (Qilin Bell, catalog number HD201805135); MiniL-12G large-capacity mini centrifuge (Maigao, catalog number MinilL-12G), etc.

[0100] (3) Primer synthesis and preparation Based on the primers synthesized in Example 2, the second round of primers used was the Novagen NM344 kit.

[0101] (4) Extraction of nucleic acid from samples The samples used in the present invention are provided by the units to which this application belongs, and these samples are real clinical samples that are confirmed to contain at least one target pathogen. Add internal standard to the cerebrospinal fluid sample, mix, centrifuge, and discard the supernatant (about 400 μL of residual liquid). After the residual liquid is broken, lysate LE and proteinase K are added thereto, and after mixing by blowing, about 600 μL of supernatant is drawn after brief centrifugation and added to the 1 / 7 row of holes in a 96-well deep well plate, and nucleic acid extraction is completed using a self-developed nucleic acid extraction kit. The concentration of the extracted nucleic acid is then determined using the Qubit dsDNA reagent.

[0102] (5) PCR amplification Reverse transcription kit (Abclonal), one-round amplification kit (Takara) and two-round amplification kit (Vazyme) were used to synthesize cDNA and prepare libraries for the nucleic acids involved above. Qubit dsDNA reagent was used to determine the product concentration after each round of amplification.

[0103] (6) Sequencing The library was prepared with nanoballs (DNB) and sequenced using the DNBSEQ-E25RS high-throughput sequencing reagent set. The sequencing was single-end 50bp sequencing. The libraries were mixed at equal mass, and DNB was prepared with 100 fmol input according to the instrument instructions. Then 20 ul of DNB termination buffer was added immediately, and the flared pipette tip was slowly pipetted to mix. The single-stranded DNA concentration of DNB was detected using the QubitssDNA detection kit. After meeting 4 ng / μL ~ 40 ng / μL, E25 sequencing was performed.

[0104] (7) Data analysis After the machine is opened, the raw data is a FastQ file, which is analyzed by an automated bioinformatics process. The FastQ file generates results after quality control, primer matching, and Kraken 2 species identification process. The specific counting method for a target pathogen-specific Reads is: all Reads match the upstream primers and UMI sequences of the target pathogen respectively, and identify the species in Kraken 2. Reads that are 100% matched with the upstream primers of the target pathogen and correctly report the target pathogen or target pathogen atomic level in Kraken 2 are then corrected by the number of UMI species connected to it and the number of cleanreads after quality control filtering of the data set. If the corrected number of Reads of a species is >= 100 times the number of Reads of the species in the blank sample, the species is considered to be positive in the sample. The calculation formula for the corrected number of Reads is as follows: (UMI-normalizedReads / cleanreads)*1000000 (8) Correction Reads Results of Target Pathogens: The qualified primers of the three target pathogens were mixed into one system, and the above steps were followed to test the clinical samples containing one of the target pathogens. A total of three positive samples (samples 1-3) and one blank control sample (sample 4) were tested. The test results are shown in Table 16. The positive samples only show the number of corrected reads of the species known to be contained therein, and the blank control shows the number of corrected reads of the three target pathogens for the convenience of calculation: Table 16 Detection results of the target pathogen primer mixture system in real clinical samples

[0105] In summary, the hybrid system detected all target species contained in three positive real clinical samples, verifying the effectiveness of the primers designed and optimized by the present invention in actual clinical applications.

Claims

1. A primer design method for multiplex PCR targeted sequencing technology, characterized in that: The method comprises the following steps in sequence: S100: Determine the best reference genome and construct a high-quality genome library of the target pathogen; S200: constructing a phylogenetic tree for the high-quality genome library; S300: traversing the optimal reference genome with a fixed step length to obtain candidate primer pairs; S400: Evaluate the candidate primer pairs, including scoring of coverage, specificity and dimer binding.

2. The primer design method for multiplex PCR targeted sequencing technology according to claim 1, characterized in that: The step S100 specifically includes the following steps: S110: Obtain the best reference genome of the target pathogen; S120: obtaining a number of genomes of other species with the most recent sequencing time, wherein the genomes of other species are obtained in order of priority from high to low assembly levels; S130: Screening out other genomes of the same species that have high completeness and meet the sample quantity requirements. When the sample quantity is insufficient, preferentially extending the sequencing time of step S120.

3. The primer design method for multiplex PCR targeted sequencing technology according to claim 2, characterized in that: The method for determining the best reference genome is: selecting a best reference genome from the RefSeq database or the GenBank database in the priority order of Representative Genome > Reference Genome > Complete Genome > Chromosome > Scaffold > Contig; When there are multiple genomes or genome structures at the same level, the genome or genome structure with the most recent sequencing year and the highest CheckM value is selected; Among them, the priority of sequencing year is greater than CheckM value; and if the sequencing years and CheckM values ​​of multiple genomes or genome structures of the same level are the same, one is selected arbitrarily.

4. The primer design method for multiplex PCR targeted sequencing technology according to claim 2 or 3, characterized in that: For pathogens that are taxonomically classified as "species", the number of samples shall not be less than 10; For pathogens taxonomically classified as "species complex" or "genus", the sample size is 50.

5. The primer design method for multiplex PCR targeted sequencing technology according to claim 4, characterized in that: CheckM2 was used to perform integrity testing and select genomes with integrity greater than or equal to 95%.

6. The primer design method for multiplex PCR targeted sequencing technology according to claim 1, characterized in that: The genome evolution tree obtained in step S200 meets the following requirements: the branches of the genome evolution tree are uniform, there are no isolated branches with obviously long evolutionary distances, and there is no situation where most of the genomes are densely distributed in the same branch.

7. The primer design method for multiplex PCR targeted sequencing technology according to claim 1, characterized in that: In step S300, the whole genome is sliced ​​in a shingled manner with a sliding window length of 300 bp and a sliding step length of 20 bp, and primers are designed on each fragment using primer3 software.

8. A primer design method for multiplex PCR targeted sequencing technology according to any one of claims 1-3, 5-7, characterized in that: The coverage evaluation method is to match the candidate primer pairs to each other genome of the same species in the high-quality genome library and calculate the primer coverage: Primer coverage = number of qualified matching genomes / total number of genomes × 100%.

9. The primer design method for multiplex PCR targeted sequencing technology according to claim 8, characterized in that: A qualified match satisfies at least the following conditions simultaneously: both the upstream primer and the downstream primer can achieve a match, the matching length of the primer amplicon is not less than 95% of the total length of the amplicon, and the primer amplicon length is 100 - 300 bp.

10. A primer design method for multiplex PCR targeted sequencing technology according to any one of claims 1-3, 5-7 or 9, characterized in that: The specific evaluation method comprises the following steps: S421: aligning the candidate primer pairs in pairs in the NCBI nt library against all genomes in the library; S422: After alignment, a primer insert fragment is obtained, wherein the primer insert fragment refers to a nucleic acid sequence defined and amplified by a pair of complementary oligonucleotide primers during PCR, wherein the starting end of the primer insert fragment is defined by the 5' end of the forward primer, and the ending end is defined by the 3' end of the reverse primer; S423: inputting the primer insert fragment into Kraken 2 software for verification, thereby screening primer pairs with qualified specificity; S424: When the species of the genome to which the primer insert fragment belongs is consistent with the primer target species, the species discrimination result of Kraken 2 must be that species; or when the species of the genome to which the primer insert fragment belongs is inconsistent with the primer target species, the species discrimination result of Kraken 2 must be consistent with the genome species.

11. The primer design method for multiplex PCR targeted sequencing technology according to claim 10, characterized in that: The dimer binding score is: ； When any one of the candidate primer pairs No less than 10 3 When , the primer pair corresponding to the primer is screened out; Wherein, GC-count indicates the number of GC bases in the dimer binding region, L indicates the base length of the dimer binding region; d1 indicates the distance from the starting point of the dimer binding region to the 3' end of the nearest primer, and d2 indicates the distance from the end point of the dimer binding region to the 3' end of the nearest primer.

12. A primer design method for multiplex PCR targeted sequencing technology according to claim 1-3, 5-7, 9 or 11, characterized in that: The method is used for multiplex PCR targeted pathogen detection.

13. The primer design method for multiplex PCR targeted sequencing technology according to claim 12, characterized in that: The pathogens include bacteria, fungi and viruses.

Images