STR typing method applied to next-generation sequencing data

Abstract

The invention provides an STR typing method applied to next-generation sequencing data, and relates to the technical field of biology. The method specifically comprises the following steps: customizing an STR locus configuration file; constructing a flanking sequence k-mer index database of the STR gene locus; an STR allele sequence is extracted from each READ of the sequence comparison file in combination with the CIGAR value and the k-mer index; constructing a PCR shadow peak model of the STR gene locus, and estimating parameters of the PCR shadow peak model by using a maximum likelihood estimation method; estimating a candidate allele frequency by using the next-generation sequencing sample data set; and calculating posterior probabilities of all candidate alleles by using an allele set extracted from a given sequence comparison file, and deducing the most probable allele. The STR typing method has the advantages of higher detection rate, higher typing accuracy and higher typing speed, and can be used for whole genome sequencing data and high-coverage targeted sequencing data or amplicon sequencing data.

Description

Technical Field

[0001] The present invention relates to the field of biotechnology, and in particular to a STR typing method applied to second-generation sequencing data. Background Art

[0002] STR (Short Tandem Repeats), also known as microsatellites, is a sequence of up to 100 nucleotides in length, consisting of a series of 2-6 nucleotide DNA repeating units. Studies have shown that, unlike SNVs (Single Nucleotide Variations) which have only two alleles, STRs are highly polymorphic in the human genome, and are therefore widely used in medical genetics, population genetics, and forensic science. In the medical field, STRs have been shown to play a key role in more than 65 Mendelian genetic diseases, and are involved in gene regulation and the formation of complex traits. In forensic science, STRs are widely used in individual identification, paternity testing, and surname inference.

[0003] STR allele typing is mainly based on capillary electrophoresis (CE) and fluorescent labeling multiple amplification STR typing methods. This method is simple to operate, economical, and fast to detect. Therefore, it has been used as the "gold standard" in forensic DNA laboratories for decades. However, STR typing technology based on capillary electrophoresis is limited by spectral resolution, resulting in only typing dozens of STR loci at a time, which cannot meet the needs of large-scale samples for high-throughput STR typing.

[0004] With the emergence and development of high-throughput sequencing technology, researchers can simultaneously amplify hundreds of STR loci through multiplex amplicon sequencing (MAS), which can largely eliminate the limitation of the number of STR loci typing by traditional CE technology. In addition, STR typing technology based on high-throughput sequencing technology has higher sensitivity and accuracy. Studies have shown that 62pg of genomic DNA can successfully type 97.7% of alleles. However, due to multiple factors such as sequencing depth, flanking sequence repeatability, and PCR shadow peaks, it is difficult to obtain accurate STR typing results from low-coverage sequencing data. At present, there are many STR typing methods based on high-throughput sequencing data, such as STRait Razor, HipSTR, STRling, ExpansionHunter, STRinNGS, etc., but these methods have problems such as poor compatibility, low performance, and low accuracy when performing STR typing. STRait Razor and HitSTR cannot handle STR loci with highly repetitive flanking sequences well; STRling and ExpansionHunter are designed for long-fragment STR loci and cannot work well in single-end sequencing data; STRinNGS matches STR repeat units through regular expressions, resulting in slow typing speed and low typing accuracy.

[0005] For example, Chinese patent CN113362892A discloses a method for detecting and typing the number of short tandem repeat sequences. This method relies heavily on the conservatism of the flanking sequences. If there are variations in the flanking sequences, it will have a significant impact on the comparison, causing problems with STR boundary identification, further affecting STR typing. Chinese patent CN117037906A discloses a method for typing short tandem repeat sequences based on second-generation sequencing. This method relies heavily on the reference genome. If the test sample differs greatly from the reference genome, the reads will be difficult to match to the reference genome, and even if they are matched, their boundaries will not be ready to match.

[0006] Based on the above problems, developing a STR typing method for second-generation sequencing data to improve its typing efficiency and accuracy and provide basic support for the accurate typing of massive omics data is the research focus of researchers in this field. Summary of the invention

[0007] In view of the above problems, the present invention provides a STR typing method applied to second-generation sequencing data, which can obtain the alleles of a specified STR locus and their likelihood probabilities by typing from sequencing data.

[0008] To achieve the above purpose, the technical solution adopted by the present invention is as follows:

[0009] In one aspect, the present invention provides a STR typing method applied to second-generation sequencing data, comprising the following steps:

[0010] S1: Customized STR locus profile;

[0011] S2: Construct k-mer index of flanking sequences of STR loci: construct a k-mer index library for the flanking sequences of each STR locus in the configuration file;

[0012] S3: for each STR locus in the STR locus configuration file of step S2, execute steps S4-S7 respectively;

[0013] S4: Extract STR locus allele sequence: extract STR allele sequence from each READ of the sequence alignment file by combining CIGAR value and k-mer index to generate the allele set of STR locus;

[0014] S5: Construct a PCR shadow peak model for the STR locus, and use the next-generation sequencing sample data set and the maximum likelihood estimation method to estimate the parameters of the PCR shadow peak model;

[0015] S6: Estimation of candidate allele frequencies at STR loci: Estimate candidate allele frequencies using a next-generation sequencing sample dataset;

[0016] S7: Infer the most likely allele and its posterior probability of the STR locus in the sequence alignment file: Use the allele set extracted in S4 to calculate the posterior probability of all candidate alleles to infer the most likely allele.

[0017] Preferably, in S1, the STR locus configuration file is used to define the main information of the STR locus and the precise position on the reference genome.

[0018] Preferably, in S1, the STR locus configuration file comprises the following fields: STR locus name, chromosome, start position, end position, repeat unit length, number of excluded bases and whether it is haplotype.

[0019] Preferably, the repeat unit length is used to calculate allele size.

[0020] Preferably, the excluded base numbers can exclude sequences that are not required to participate in the calculation of alleles.

[0021] Preferably, the haplotype is used to genotype a STR locus consisting of two segments in a sex chromosome.

[0022] Preferably, in S2, the flanking sequence is a base sequence connected to both sides of the 5' and 3' ends of the STR core sequence.

[0023] Preferably, in S2, the k-mer index is a local unique k-mer subsequence, which can improve the matching accuracy of the STR locus flanking sequence.

[0024] Preferably, in S2, the construction of the STR locus flanking sequence k-mer index adopts a dynamic k-mer index construction strategy, that is, according to the repeatability of the STR locus flanking sequence, the optimal length k is selected as the k-mer subsequence length of the locus, and then the k-mer subsequence length k and the relative position p of the k-mer in the flanking sequence are stored in a hash table to prepare for the subsequent extraction of alleles.

[0025] Preferably, in S4, the alignment file is a BAM file after alignment to the reference genome.

[0026] Preferably, in S4, the specific operation of extracting the STR allele sequence is:

[0027] S4-1, using the sequence alignment file as input, first calculate the STR flanking sequence position using the CIGAR value of each READ and the position of the corresponding STR locus on the reference genome;

[0028] S4-2, then use the flanking sequence in the configuration file to verify the calculated READ flanking sequence, so as to directly calculate the STR allele;

[0029] S4-3. When the insertion or deletion in the CIGAR value appears at the junction of the flanking sequence and the core sequence and the number of bases does not exceed 2 times the length of the STR repeat unit, the CIGAR value correction method is used to move the insertion or deletion in the 5' and 3' flanking sequences that are identical to the repeat unit and have an incorrect alignment back to the STR core region, and finally the STR allele sequence is extracted using the corrected CIGAR value;

[0030] S4-4. When the allele length of the sample on the STR genome exceeds 40% of the allele length of the reference genome at that site, the sequence alignment software often introduces completely wrong CIGAR values. In order to make up for the defects of the CIGAR method, the k-mer index method is enabled to extract alleles: for each READ sequence, the corresponding STR site information is extracted from the k-mer index library constructed in S2, and then the flanking sequences are matched and the STR alleles are calculated.

[0031] Preferably, in S5, the PCR shadow peak model comprises parameters u and d, wherein u represents the probability of having one more repeat unit compared to the true allele, and d represents the probability of having one less repeat unit compared to the true allele.

[0032] Preferably, in S5, the distribution of the PCR shadow peak model is as follows:

[0033] For a given STR locus, x represents the number of repeat units that differ between the observed allele and the true allele, which can be divided into three cases:

[0034] Same as the true allele (x = 0);

[0035] One more repeat unit compared to the true allele (x = 1);

[0036] and one less repeat unit compared to the true allele (x = -1);

[0037] Therefore, the distribution of the PCR shadow peak model can be written as:

[0038]

[0039] Preferably, in S5, the specific operation of estimating the parameters of the PCR shadow peak model using the second-generation sequencing sample data set and the maximum likelihood estimation method is:

[0040] S5-1. For a given individual i, there are n i If there are reads that completely span the locus, the PCR shadow peak set is in represents the difference between the observed allele and the true allele in the j-th READ across the current locus for individual i;

[0041] Then, the likelihood function for N individuals in the second-generation sequencing sample data set is as follows:

[0042]

[0043] S5-2. Based on formulas (1) and (2), the likelihood function of u and d can be calculated:

[0044]

[0045] in:

[0046] T represents the sum of the number of reads of N samples that completely span the current STR locus;

[0047] m1 indicates the number of alleles in the reads that have one more repeat unit than the true allele;

[0048] m2 indicates the number of alleles in the reads that have one less repeat unit compared to the true allele;

[0049] S5-3. By using formulas (3) and (4), we can estimate the maximum likelihood estimates of parameters u and d:

[0050]

[0051] Preferably, in S5, for estimating the user-defined PCR shadow peak model parameters, the present invention needs to screen the input samples through strict standards to estimate and obtain reliable PCR shadow peak model parameters, and its basic rules include:

[0052] (1) The number of reads that completely span the STR locus must exceed the set minimum limit (the default is 10);

[0053] (2) For the allele with the largest number of reads support, its read support number must exceed 70% of the total number of reads at the STR locus.

[0054] Preferably, in S6, the estimating of the candidate allele frequency specifically comprises: calculating the current STR locus allele frequency distribution through the second generation sequencing sample data set input by the user, and recalculating the candidate allele frequency of the current locus when the sample size increases.

[0055] Preferably, in S6, the standard for estimating the candidate allele frequency is:

[0056] (1) The number of sequencing samples used to estimate the allele frequency distribution is no less than 20;

[0057] (2) The candidate allele must appear in at least two samples.

[0058] Preferably, in S6, for a data set where the sequencing samples are insufficient to estimate the allele frequencies of the custom STR loci, the same frequency value is assigned to each potential STR allele that occurs, thereby eliminating the influence of the prior probability on the posterior probability estimate.

[0059] Preferably, in said S7, in S7, the most likely allele of the current STR locus in the sequence alignment file and its posterior probability are inferred by:

[0060] S7-1. If the current STR locus is defined as a haplotype in the configuration file, and there are n reads that completely span the current STR locus, the allele set composed of the alleles observed in these reads can be expressed as A = {a1, a2, ..., an}, then, the posterior probability of each possible candidate allele a(a∈A) can be calculated by the following formula:

[0061]

[0062] Among them, f a represents the frequency of candidate allele a in the current STR locus;

[0063] a i represents the i-th observed allele of the sample at the current locus;

[0064] a+1 and a-1 represent the shadow peak alleles with one more repeat unit and one less repeat unit compared with the given candidate allele a, respectively;

[0065] others indicates the observed allele a i A number of more or fewer repeat units than a given candidate allele a;

[0066] S7-2. If the current STR locus is defined as non-haploid in the configuration file. In order to estimate the posterior probability of two candidate allele combinations, assume that the two alleles of the STR locus are j (j∈A) and k (k∈A); then the joint posterior probability of all pairwise combinations of alleles can be estimated:

[0067] P(allele=(j,k)|A,u,d)∝#

[0068]

[0069] Among them, f j ,f k represents the frequency of candidate allele j and allele k in the current STR locus;

[0070] a i represents the i-th observed allele of the sample at the current locus;

[0071] a+1 and a-1 represent the shadow peak alleles with one more repeat unit and one less repeat unit compared with the given candidate allele a, respectively;

[0072] S7-3. Obtain the posterior probabilities of all candidate alleles or allele combinations of the current STR locus through calculation, and select the allele or allele combination with the largest posterior probability as the most likely allele or allele combination at the current STR locus in the sequence alignment file.

[0073] Compared with the prior art, the present invention has the following beneficial effects:

[0074] The present invention proposes a STR typing technology solution that can be applied to second-generation sequencing data. It can type the alleles of a specified STR locus and their likelihood probabilities from the sequencing data, and has higher detection rate and typing accuracy. In addition to being applicable to whole genome sequencing (WGS) sequencing data, it can also be applied to high-coverage targeted sequencing data or amplicon sequencing data. The STR typing speed is greatly improved through algorithm optimization. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 A framework diagram of the STR typing method applied to second-generation sequencing data provided by the present invention.

[0076] Figure 2 A technical roadmap for the STR typing method applied to second-generation sequencing data provided by the present invention. DETAILED DESCRIPTION

[0077] In order to make the technical means, creative features, purpose and effect of the present invention easy to understand, the present invention is further explained below in conjunction with specific embodiments, but the following embodiments are only preferred embodiments of the present invention, not all. Based on the embodiments in the implementation mode, other embodiments obtained by those skilled in the art without making creative work all belong to the protection scope of the present invention. It is worth noting that the raw materials used in the present invention are all common commercial products, and their sources are not specifically limited. The technology and scientific terms used in the embodiments have the meanings commonly understood by those of ordinary skill in the art to which the present invention belongs.

[0078] Example 1

[0079] A STR typing method applied to next-generation sequencing data comprises the following steps:

[0080] (1) Step 1: Customize the STR locus profile to determine the necessary parameters such as the position of each STR locus on the reference genome, the length of the repeat unit, and the number of excluded bases

[0081] The STR configuration file is used to define the main information of the STR locus and its precise position on the reference genome. It mainly contains 7 fields, namely, STR locus name, chromosome, start position, end position, repeat unit length, number of excluded bases, and whether it is haplotype. The length of the repeat unit is used to calculate the allele size; the number of excluded bases can exclude sequences that do not need to participate in the allele calculation; whether it is haplotype is used to type the STR loci composed of two fragments in the sex chromosome (such as DYS385ab).

[0082] (2) Step 2: Construct a k-mer index of the flanking sequences of the STR loci. Construct a k-mer index for the flanking sequences (5' and 3') of each STR locus in the configuration file for subsequent positioning of the flanking sequences.

[0083] The extraction of STR alleles must accurately match the 5' and 3' flanking sequences of the STR core region, so this method introduces the k-mer method to match the flanking sequences. Since the similarity of the flanking sequences of many STR loci is high, if the k-mer index is directly constructed for the flanking sequence, the subsequent STR locus matching error rate will be greatly increased. For this reason, the method uses a local unique k-mer subsequence as the flanking sequence index to improve the STR flanking sequence matching accuracy. In addition, due to the large difference in the repetitiveness of the flanking sequences of different STR loci, the use of a fixed-length k-mer index cannot cope with this situation well. For this reason, the present invention adopts a dynamic k-mer index construction strategy, that is, the optimal length k is selected as the k-mer subsequence length of the locus according to the repetitiveness of the STR locus flanking sequence. Finally, the k-mer subsequence length k and the relative position p of the k-mer in the flanking sequence are stored in a hash table to prepare for the subsequent allele extraction.

[0084] (3) Step 3: Extract STR allele sequences. Extract STR allele sequences from each READ in the sequence alignment file (BAM) by combining CIGAR value and k-mer index.

[0085] The extraction of STR allele sequence needs to locate the STR flanking sequence first. The present invention combines CIGAR value and k-mer index to locate the flanking sequence in READ, so as to improve the success rate of STR allele sequence matching. This method takes the BAM file after alignment to the reference genome as input, first uses the CIGAR value of each READ and the position of the corresponding STR locus on the reference genome to calculate the STR flanking sequence position, and then uses the flanking sequence in the configuration file to verify the calculated READ flanking sequence, so as to directly calculate the STR allele size. However, when the flanking sequence is highly similar to the core sequence, it will cause the sequence alignment tool to indicate the position of the insertion or deletion with an incorrect CIGAR value. If allele typing is performed directly by this CIGAR value, it will often lead to incorrect typing results. In order to solve this problem, the present invention adopts the CIGAR value correction method to move the insertion or deletion identical to the repeat unit in the 5' and 3' end flanking sequences back to the STR core region, and finally uses the corrected CIGAR value to extract the STR allele sequence.

[0086] When the allele length of a sample on the STR genome exceeds 40% of the allele length of the reference genome at that site, the sequence alignment tool will mark a portion of the sample READ sequence as softclip to obtain a higher alignment score. However, this alignment strategy will lead to the introduction of completely wrong CIGAR values, resulting in errors in STR typing results. In order to make up for the defects of the CIGAR method, the present invention introduces a flanking sequence matching strategy based on the k-mer method, that is, when the CIGAR method fails, the k-mer method is enabled to extract alleles. For each READ sequence, the corresponding STR site information (including the k-mer length k, the position of the k-mer in the flanking sequence, etc.) is extracted from the k-mer index library constructed in step (2), and then the flanking sequences are matched and the STR allele size is calculated.

[0087] (4) Step 4: Construct a PCR shadow peak model for the STR locus and use the sequencing sample data provided by the user to estimate the probability that the PCR shadow peak has one more or one less repeat unit than the true allele.

[0088] PCR shadow peaks have a great impact on the accurate typing of STR, especially in low-coverage sequencing data. The PCR shadow peak model can significantly improve the accuracy of STR allele typing. Since the distribution of PCR shadow peaks of different STR loci varies greatly, the present invention introduces a corresponding PCR shadow peak distribution for each STR locus to eliminate the impact of PCR shadow peaks on STR typing through maximum a posteriori probability estimation. The PCR shadow peak model contains two parameters u and d, where u represents the probability of having one more repeat unit than the true allele, and d represents the probability of having one less repeat unit than the true allele. Therefore, for a given STR locus, x represents the number of repeat units that differ between the observed allele and the true allele, which is divided into three cases, namely, the same size as the true allele (x=0), one more repeat unit than the true allele (x=1), and one less repeat unit than the true allele (x=-1). Therefore, the PCR shadow peak distribution can be written as:

[0089]

[0090] After obtaining the PCR shadow peak model, this method needs to estimate the model parameters through the sequencing reads of a given sample. i If there are reads that completely span the locus, the PCR shadow peak set is in It represents the difference between the observed allele and the true allele in the jth READ across the current locus of individual i. Then, the likelihood function for n individuals in the second-generation sequencing sample data set is as follows:

[0091]

[0092] Based on formulas (1) and (2), the maximum likelihood function for u and d can be calculated:

[0093]

[0094] in:

[0095] T represents the sum of the number of reads of N samples that completely span the current STR locus;

[0096] m1 indicates the number of alleles in the reads that have one more repeat unit than the true allele;

[0097] n2 indicates the number of alleles in the reads that have one less repeat unit compared to the true allele;

[0098] S5-3. By using formulas (3) and (4), we can estimate the maximum likelihood estimates of parameters u and d:

[0099]

[0100] For the estimation of user-defined PCR shadow peaks, this method needs to screen the input samples through strict standards to estimate reliable PCR shadow peak model parameters. Its basic rules include: (1) The number of reads that completely span the STR locus must exceed the set minimum limit (the default is 10) (2) For the allele with the largest number of reads support, its read support number must exceed 70% of the total number of reads on the STR locus.

[0101] (5) Step 5: Estimate candidate allele frequencies using the provided sequencing sample data

[0102] The allele frequency of STR loci has a relatively stable distribution characteristic in the population, which plays an important role in the estimation of posterior probability during STR typing. In particular, in low-coverage sequencing data, accurate allele frequency distribution can greatly improve the accuracy of STR typing results. For some common STR loci, such as the allele frequency of forensic STR loci in the population, it can be easily obtained from the NIST (National Institute of Standards and Technology) website. However, for user-defined STR loci, it is often difficult to obtain ready-made allele frequency distribution data from the Internet. Therefore, the present invention can automatically estimate the allele frequency distribution of the customized STR loci through the sequencing sample data input by the user, and can continuously correct the locus allele frequency data as the sample size increases, thereby improving the accuracy of STR allele typing. In order to obtain accurate and reliable STR locus allele frequency distribution, the present invention sets a series of strict standards: (1) the number of sequencing samples used to estimate the allele frequency distribution is not less than 20 (2) the candidate allele must appear in at least 2 samples. For data sets where the sequencing samples are insufficient to estimate the allele frequencies of custom STR loci, the present invention assigns the same frequency value to each potential STR allele that occurs, thereby eliminating the influence of the prior probability on the posterior probability estimate.

[0103] (6) Step 6: Infer the most likely allele and its likelihood probability, and use the extracted allele set to calculate the maximum posterior probability of all candidate alleles to infer the most likely allele

[0104] Since only a small number of reads can completely span the entire STR locus in sequencing data with low coverage, coupled with the influence of PCR shadow peaks, it is difficult to accurately infer the alleles of the sample only by the number of supported reads. However, when the PCR shadow peak ratio of the STR locus and the frequency of the allele in the population are known, these prior knowledge can be used to estimate the posterior probability of the allele, thereby greatly improving the accuracy of allele typing. The present invention uses the estimated PCR shadow peak ratio of the STR locus and the STR allele frequency distribution to calculate the posterior probability of each allele combination. Finally, the allele with the largest posterior probability is used as the final STR typing result, and the maximum posterior probability estimation process is as follows:

[0105] For a sample, assuming that there are n reads that completely span the current STR locus, the allele set composed of these reads can be expressed as A = {a1, a2, ..., a nThen, the maximum posterior probability for each possible candidate allele a (a∈A) can be calculated by the following formula:

[0106] P(allele=a|A,u,d)∝#

[0107]

[0108] f in formula (6) a The frequency of candidate allele a in STR locus l; a i represents the i-th observed allele of the sample at the current locus; a+1 and a-1 represent the shadow peak alleles that increase and decrease one repeat unit compared with the given candidate allele a, respectively; others represent the observed allele a i There are several more or fewer repeat units than the given candidate allele a. Since the autosomal STR locus has two alleles, in order to estimate the posterior probability of this type of allele, this method assumes that the two alleles of the STR locus are j (j∈A) and k (k∈A); then the joint probability of all pairwise combinations of alleles can be estimated:

[0109] P(allele=(j,k)|A,u,d)∝#

[0110]

[0111] Test example:

[0112] 1. The present invention was benchmarked against three mainstream STR typing software (HipSTR, STRait Razor, and STRinNGS). The test data consisted of 3000 simulated WGS (single-end, 150-bp) sequencing data with sequencing depths of 10× (1000), 30× (1000), and 50× (1000), respectively. The test STR loci were 43 commonly used STR loci.

[0113] Among them, HipSTR refers to the literature "Willems, Thomas, Dina Zielinski, Jie Yuan, Assaf Gordon, Melissa Gymrek, and Yaniv Erlich. "Genome-wide profiling of heritableand de novo STR variations." Nature methods 14, no. 6 (2017): 590-592.";

[0114] STRait Razor reference literature "Warshauer, David H., Jonathan L. King, and Bruce Budowle. "STRait razor v2.0: the improved STR allele identification tool-razor." Forensic Science International: Genetics 14 (2015): 182-186.";

[0115] STRinNGS reference literature" Carina Xiaoqin Qian, Halimureti, Sima Yijiang, and Claus "STRinNGS v2.0: Improved tool for analysis and reporting of STR sequencing data." Forensic Science International: Genetics 48(2020):102331."

[0116] The test results show that the typing method of the present invention has a 100% detection rate and a 99.37% typing accuracy rate in common WGS (30×) simulated sequencing data. The accuracy rates of the other three methods mentioned above are 94.85%, 86.90% and 77.28% respectively.

[0117] 2. In addition to being applicable to low-coverage WGS sequencing data, the present invention can also be applied to high-coverage targeted sequencing data or amplicon sequencing data. In order to verify its typing effect, the present invention uses 687 real targeted sequencing data (single-end, 400-bp), the average sequencing depth is 1570×, and the test STR loci are 53 commonly used STR loci.

[0118] The test results show that the detection rate and accuracy of the typing method of the present invention are as high as 99.64% and 99.9%. In comparison, the detection rates and accuracy of the other three methods mentioned above are 80.26% / 85.31%, 98.51% / 99.65%, and 97.32% / 96.98, respectively.

[0119] 3. The present invention greatly improves the STR typing speed through algorithm optimization. Using 10 high-coverage real sequencing data (average 3216×) for typing speed testing, the present invention only needs 136 seconds to complete STR allele typing, and its running speed is 79 times that of HipSTR and more than 10,000 times that of STRinNGS.

[0120] In summary, the STR typing method specifically applied to second-generation sequencing data of the present invention has a higher detection rate and typing accuracy. In addition to being applicable to low-coverage WGS sequencing data, it can also be applied to high-coverage targeted sequencing data or amplicon sequencing data, and the STR typing speed is faster.

[0121] Finally, it should be noted that the above content is only used to illustrate the technical solution of the present invention, rather than to limit the scope of protection of the present invention. Simple modifications or equivalent substitutions of the technical solution of the present invention by ordinary technicians in this field do not deviate from the essence and scope of the technical solution of the present invention.

Claims

1. A STR typing method applied to second-generation sequencing data, characterized in that: The following steps are involved: S1: Customized STR locus profile; S2: Construct k-mer index of flanking sequences of STR loci: construct a k-mer index library for the flanking sequences of each STR locus in the configuration file; S3: for each STR locus in the STR locus configuration file of step S2, execute steps S4-S7 respectively; S4: Extract STR locus allele sequence: extract STR allele sequence from each READ of the sequence alignment file by combining CIGAR value and k-mer index to generate the allele set of STR locus; S5: Construct a PCR shadow peak model for the STR locus, and use the next-generation sequencing sample data set and the maximum likelihood estimation method to estimate the parameters of the PCR shadow peak model; S6: Estimation of candidate allele frequencies at STR loci: Estimate candidate allele frequencies using a next-generation sequencing sample dataset; S7: Infer the most likely allele and its posterior probability of the STR locus in the sequence alignment file: Use the allele set extracted in S4 to calculate the posterior probabilities of all candidate alleles to infer the most likely allele.

2. The STR typing method according to claim 1, characterized in that In S1, the STR locus configuration file includes the following fields: STR locus name, chromosome, start position, end position, repeat unit length, number of excluded bases, and whether it is haplotype.

3. The STR typing method according to claim 1, characterized in that: In S2, the k-mer index is a locally unique k-mer subsequence.

4. The STR typing method according to claim 1, characterized in that: In S2, the k-mer index of the STR locus flanking sequence is constructed, and a dynamic k-mer index construction strategy is adopted, that is, the optimal length k is selected as the k-mer subsequence length of the locus according to the repeatability of the STR locus flanking sequence, and then the k-mer subsequence length k and the relative position p of the k-mer in the flanking sequence are stored in a hash table to prepare for the subsequent extraction of alleles.

5. The STR typing method according to claim 1, characterized in that: In S4, the specific operation of extracting the STR allele sequence is: S4-1, using the sequence alignment file as input, first calculate the STR flanking sequence position using the CIGAR value of each READ and the position of the corresponding STR locus on the reference genome; S4-2, then use the flanking sequence in the configuration file to verify the calculated READ flanking sequence, so as to directly calculate the STR allele; S4-3. When the insertion or deletion in the CIGAR value appears at the junction of the flanking sequence and the core sequence and the number of bases does not exceed 2 times the length of the STR repeat unit, the CIGAR value correction method is used to move the insertion or deletion in the 5' and 3' flanking sequences that are identical to the repeat unit and have an incorrect alignment back to the STR core region, and finally the STR allele sequence is extracted using the corrected CIGAR value; S4-4. When the allele length of the sample on the STR genome exceeds 40% of the allele length of the reference genome at that site, the sequence alignment software often introduces completely wrong CIGAR values. In order to make up for the defects of the CIGAR method, the k-mer index method is enabled to extract alleles: for each READ sequence, the corresponding STR site information is extracted from the k-mer index library constructed in S2, and then the flanking sequences are matched and the STR alleles are calculated.

6. The STR typing method according to claim 1, characterized in that: In S5, the PCR shadow peak model includes parameters u and d, where u represents the probability of having one more repeat unit compared to the true allele, and d represents the probability of having one less repeat unit compared to the true allele.

7. The STR typing method according to claim 1, characterized in that: In S5, the distribution of the PCR shadow peak model is as follows: For a given STR locus, x represents the number of repeat units that differ between the observed allele and the true allele, which can be divided into three cases: Same as the true allele (x = 0); One more repeat unit compared to the true allele (x = 1); and one less repeat unit compared to the true allele (x = -1); Therefore, the distribution of the PCR shadow peak model can be written as:

8. The STR typing method according to claim 1, characterized in that: In S5, the specific operation of estimating the parameters of the PCR shadow peak model using the second-generation sequencing sample data set and the maximum likelihood estimation method is: S5-1. For a given individual i, there are n i If there are reads that completely span the locus, the PCR shadow peak set is in represents the difference between the observed allele and the true allele in the j-th READ across the current locus for individual i; Then, the likelihood function for N individuals in the second-generation sequencing sample data set is as follows: S5-2. Based on formulas (1) and (2), the likelihood function of u and d can be calculated: in: T represents the sum of the number of reads of N samples that completely span the current STR locus; m1 indicates the number of alleles in the reads that have one more repeat unit than the true allele; m2 indicates the number of alleles in the reads that have one less repeat unit compared to the true allele; S5-3. By using formulas (3) and (4), we can estimate the maximum likelihood estimates of parameters u and d:

9. The STR typing method according to claim 1, characterized in that: In S6, the estimation of the candidate allele frequency is specifically as follows: statistically analyzing the current STR locus allele frequency distribution through the second generation sequencing sample data set input by the user; and re-stating the candidate allele frequency of the current locus when the sample size increases.

10. The STR typing method according to claim 1, characterized in that: In S7, the most likely allele and its posterior probability of the current STR locus in the sequence alignment file are inferred. The specific operations are: S7-1. If the current STR locus is defined as a haplotype in the configuration file, and there are n reads that completely span the current STR locus, the allele set consisting of the alleles observed in these reads can be expressed as A = {a1, a2, …, a n }, then, the posterior probability of each possible candidate allele a(a∈A) can be calculated by the following formula: Among them, f a represents the frequency of candidate allele a in the current STR locus; a i represents the i-th observed allele of the sample at the current locus; a+1 and a-1 represent the shadow peak alleles with one more repeat unit and one less repeat unit compared with the given candidate allele a, respectively; others indicates the observed allele a i A number of more or fewer repeat units than a given candidate allele a; S7-2. If the current STR locus is defined as non-haplotype in the configuration file, in order to estimate the posterior probability of two candidate allele combinations, assume that the two alleles of the STR locus are j (j∈A) and k (k∈A); then the joint posterior probability of all pairwise combinations of alleles can be estimated: P(allele=(j,k)|A,u,d)∝# Among them, f j ,f k represents the frequency of candidate allele j and allele k in the current STR locus; a i represents the i-th observed allele of the sample at the current locus; a+1 and a-1 represent the shadow peak alleles with one more repeat unit and one less repeat unit compared with the given candidate allele a, respectively; S7-3. Obtain the posterior probabilities of all candidate alleles or allele combinations of the current STR locus through calculation, and select the allele or allele combination with the largest posterior probability as the most likely allele or allele combination at the current STR locus in the sequence alignment file.

Images