A method for STR typing applied to second-generation sequencing data

Abstract

The application provides a STR typing method applied to second-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 k-mer index library of a flanking sequence of the STR locus; extracting STR allele sequences from each READ of a sequence alignment file in combination with a CIGAR value and the k-mer index; constructing a PCR ghost peak model of the STR locus, and estimating parameters of the PCR ghost peak model by using a maximum likelihood estimation method; estimating candidate allele frequencies by using a second-generation sequencing sample data set; calculating posterior probabilities of all candidate alleles by using an allele set extracted from a given sequence alignment file, and inferring the most possible allele.The STR typing method has higher detection rate and typing accuracy, and has faster typing speed, and can be used for whole genome sequencing data, high-coverage targeted sequencing data or amplicon sequencing data.

Description

Technical Field

[0001] This invention relates to the field of biotechnology, and more specifically to an STR typing method for next-generation sequencing data. Background Technology

[0002] STR (Short Tandem Repeats), also known as microsatellite sequences, are sequences up to 100 nucleotides in length, composed of a series of DNA repeating units of 2-6 nucleotides. Studies have shown that, unlike SNVs (Single Nucleotide Variations), which have only two genotypes, STRs exhibit high polymorphism in the human genome, and are therefore widely used in medical genetics, population genetics, and forensic science. In medicine, STRs have been shown to play a crucial role in over 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 deduction.

[0003] STR allele typing is primarily based on capillary electrophoresis (CE) and fluorescently labeled multiplex amplification (FLAA) STR typing methods. These methods are simple to operate, cost-effective, and fast, and have therefore remained the "gold standard" for forensic DNA laboratories for decades. However, capillary electrophoresis-based STR typing technology is limited by spectral resolution, resulting in the typing of only a few dozen STR loci per run, which cannot meet the high throughput requirements of large-scale sample typing.

[0004] With the emergence and development of high-throughput sequencing technology, researchers can simultaneously amplify hundreds or even thousands of STR loci using multiplex amplicon sequencing (MAS), largely eliminating the limitation of traditional CE technology on the number of STR loci that can be genotyped. Furthermore, STR genotyping techniques based on high-throughput sequencing offer higher sensitivity and accuracy; studies have shown that 97.7% of alleles can be successfully genotyped from just 62 pg of genomic DNA. However, due to factors such as sequencing depth, flanking sequence repetition, and PCR shadow peaks, obtaining accurate STR genotyping results from low-coverage sequencing data is quite challenging. Currently, many STR genotyping methods based on high-throughput sequencing data exist, such as STait Razor, HipSTR, STRling, ExpansionHunter, and STRINGS; however, these methods suffer from poor compatibility, low performance, and low accuracy in STR genotyping. STRait Razor and HitSTR cannot handle STR loci with high repetition in flanking sequences well; STRling and ExpansionHunter are designed for long STR sites and cannot work well in single-end sequencing data; STRinNGS matches STR repeat units using regular expressions, resulting in slow typing speed and low typing accuracy.

[0005] For example, Chinese patent CN113362892A discloses a method for detecting and genotyping short tandem repeat sequences. This method relies heavily on the conservation of flanking sequences. If there are variations in the flanking sequences, it will greatly affect alignment, causing problems with STR boundary identification and further affecting STR genotyping. Chinese patent CN117037906A discloses a genotyping method for short tandem repeat sequences based on next-generation sequencing. This method relies heavily on the reference genome. If the test sample differs significantly from the reference genome, the reads will be difficult to align to the reference genome, and even if they are aligned, their boundaries cannot be accurately matched.

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

[0007] To address the aforementioned problems, this invention provides an STR genotyping method for second-generation sequencing data, which can obtain the alleles and their likelihood probabilities for a specified STR locus from sequencing data.

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

[0009] On the one hand, this invention provides an STR typing method for second-generation sequencing data, comprising the following steps:

[0010] S1: Custom STR locus configuration file;

[0011] S2: Constructing a k-mer index for flanking sequences of STR loci: Constructing 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 described in step S2, execute steps S4-S7 respectively;

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

[0014] S5: Construct a PCR shadow peak model for STR loci, and estimate the parameters of the PCR shadow peak model using a second-generation sequencing sample dataset and the maximum likelihood estimation method.

[0015] S6: Estimating candidate allele frequencies at STR loci: Estimating candidate allele frequencies using second-generation sequencing sample datasets;

[0016] S7: Infer the most likely alleles and their posterior probabilities for STR loci in the sequence alignment file: Calculate the posterior probabilities of all candidate alleles using the allele set extracted in S4, and infer the most likely alleles accordingly.

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

[0018] Preferably, 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 a haplotype.

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

[0020] Preferably, the excluded base number can exclude sequences that do not need to be included in the allele calculation.

[0021] Preferably, whether it is haplotype is used for STR loci composed of two segments in the genotypic chromosome.

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

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

[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, the optimal length k is selected as the k-mer subsequence length of the locus based on the repetition 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.

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

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

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

[0028] S4-2. Then, the calculated READ flanking sequence is verified using the flanking sequence in the configuration file, and the STR allele is directly calculated in this way.

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

[0030] S4-4. When the allele length of a sample on the STR genome exceeds 40% of the allele length of the reference genome at that locus, sequence alignment software often introduces completely erroneous CIGAR values. To compensate for the defects of the CIGAR method, the k-mer index method is enabled for allele extraction: For each READ sequence, the corresponding STR locus 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 includes 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.

[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, and there are three possible cases:

[0034] Identical to the true allele (x=0);

[0035] It has one more repeating unit (x=1) compared to the true allele;

[0036] And it has one less repeating unit compared to the true allele (x = -1);

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

[0038]

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

[0040] S5-1. For a given individual i, n has at the current STR locus i A read that completely spans the locus will cause the PCR peaks to aggregate as follows: in This represents the difference between the observed allele and the true allele in the j-th READ that spans the current locus for individual i.

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

[0042]

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

[0044]

[0045] in:

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

[0047] m1 indicates the number of repeat units in the reads that are one more than the actual alleles;

[0048] m2 represents the number of alleles in the reads that have one less repeat unit compared to the actual alleles;

[0049] S5-3. The maximum likelihood estimates of parameters u and d can be obtained through formulas (3) and (4):

[0050]

[0051] Preferably, in S5, for estimating user-defined PCR peak model parameters, this invention requires rigorous standard screening of input samples to estimate reliable PCR peak model parameters. The basic rules include:

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

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

[0054] Preferably, in S6, estimating the candidate allele frequency specifically involves: statistically analyzing the allele frequency distribution of the current STR locus using the second-generation sequencing sample dataset input by the user. Furthermore, the candidate allele frequency of the current locus can be recalculated as the sample size increases.

[0055] Preferably, in S6, the criterion for estimating candidate allele frequencies is:

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

[0057] (2) Candidate alleles must appear in at least two samples.

[0058] Preferably, in S6, for datasets where the sequencing samples are insufficient to estimate the allele frequencies of custom STR loci, each potential STR allele will be assigned the same frequency value to eliminate the influence of prior probability on posterior probability estimation.

[0059] Preferably, in step S7, inferring the most probable allele and its posterior probability for the current STR locus in the sequence alignment file, the specific operation is as follows:

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

[0061]

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

[0063] a i This indicates the i-th observed allele at the current locus in the sample;

[0064] a+1 and a-1 represent shadow peak alleles that have one more repeat unit and one less repeat unit compared to a given candidate allele a, respectively.

[0065] "Others" represents the observed allele a. i It has a certain number of repeating units more or less than the given candidate allele a;

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

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

[0068]

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

[0070] a i This indicates the i-th observed allele at the current locus in the sample;

[0071] a+1 and a-1 represent shadow peak alleles that have one more repeat unit and one less repeat unit compared to a given candidate allele a, respectively.

[0072] S7-3. By calculation, the posterior probabilities of all candidate alleles or allele combinations at the current STR locus are obtained. The allele or allele combination with the highest posterior probability is selected as the most likely allele or allele combination in the sequence alignment file at the current STR locus.

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

[0074] This invention proposes an STR genotyping technique applicable to second-generation sequencing data, which can obtain the alleles and their likelihood probabilities of a specified STR locus from sequencing data, with higher detection and genotyping accuracy. In addition to being applicable to whole-genome sequencing (WGS) data, it can also be applied to high-coverage targeted sequencing data or amplicon sequencing data. The algorithm optimization greatly improves the speed of STR genotyping. Attached Figure Description

[0075] Figure 1 This is a framework diagram of the STR typing method for second-generation sequencing data provided by the present invention.

[0076] Figure 2 This invention provides a technical roadmap for the STR typing method applied to second-generation sequencing data. Detailed Implementation

[0077] To make the technical means, creative features, achieved objectives, and effects of this invention readily understandable, the invention is further illustrated below with specific embodiments. However, these embodiments are merely preferred embodiments and not all embodiments. Other embodiments obtained by those skilled in the art based on the embodiments described herein without creative effort are all within the scope of protection of this invention. It is worth noting that the raw materials used in this invention are all common commercially available products, and their sources are not specifically limited. The technical and scientific terms used in the embodiments have the meanings commonly understood by those skilled in the art to which this invention pertains.

[0078] Example 1

[0079] A method for STR typing of next-generation sequencing data includes the following steps:

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

[0081] The STR profile is used to define the main information of STR loci and their precise location on the reference genome. It mainly contains seven fields: STR locus name, chromosome, start position, end position, repeat unit length, number of excluded bases, and whether it is haplotype. The repeat unit length is used to calculate allele size; the number of excluded bases can exclude sequences that do not need to participate in allele calculation; whether it is haplotype is used for STR loci consisting of two segments on the genotype chromosome (such as DYS385ab).

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

[0083] STR allele extraction requires accurate matching to the 5' and 3' flanking sequences of the STR core region. Therefore, this method introduces the k-mer method for flanking sequence matching. Since many STR loci have high flanking sequence similarity, directly constructing k-mer indexes on the flanking sequences would significantly increase the subsequent STR locus matching error rate. Therefore, this method uses locally unique k-mer subsequences as flanking sequence indexes to improve the accuracy of STR flanking sequence matching. Furthermore, due to significant differences in the repetition of flanking sequences at different STR loci, a fixed-length k-mer index cannot adequately address this situation. Therefore, this invention employs a dynamic k-mer index construction strategy, selecting the optimal length k as the k-mer subsequence length for that locus based on the repetition of the STR locus flanking sequences. Finally, the k-mer subsequence length k and the relative position p of the k-mer in the flanking sequences are stored in a hash table to prepare for 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] Extraction of STR allele sequences requires first locating the STR flanking sequences. This invention combines CIGAR values ​​and k-mer indexing to locate flanking sequences in READ sequences, thereby improving the success rate of STR allele sequence matching. This method uses the BAM file aligned to the reference genome as input. First, it calculates the STR flanking sequence positions using the CIGAR value of each READ sequence and the position of the corresponding STR locus on the reference genome. Then, it verifies the calculated READ flanking sequences using the flanking sequences in the configuration file, thus directly calculating the STR allele size. However, when the flanking sequences have high similarity to the core sequence, sequence alignment tools may use incorrect CIGAR values ​​to represent insertion or deletion positions. Directly using these CIGAR values ​​for allele typing often leads to incorrect typing results. To solve this problem, this invention employs a CIGAR value correction method to move insertions or deletions identical to repeat units in the 5' and 3' flanking sequences back to the STR core region, and finally uses the corrected CIGAR values ​​to extract the STR allele sequences.

[0086] When the allele length of a sample on the STR genome exceeds 40% of the allele length of the reference genome at that locus, the sequence alignment tool will mark a portion of the READ sequence of the sample as a softclip to obtain a higher alignment score. However, this alignment strategy can lead to the introduction of completely erroneous CIGAR values, resulting in incorrect STR genotyping results. To compensate for the shortcomings of the CIGAR method, this 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 used to extract alleles. For each READ sequence, the corresponding STR locus information (including 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 to calculate the STR allele size.

[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 actual allele.

[0088] PCR shadow peaks have a significant impact on the accurate genotyping of STRs, especially in low-coverage sequencing data. PCR shadow peak models can significantly improve the accuracy of STR allele genotyping. Since the PCR shadow peak distribution varies considerably among different STR loci, this invention introduces a corresponding PCR shadow peak distribution for each STR locus to eliminate the influence of PCR shadow peaks on STR genotyping through maximum a posteriori probability estimation. This PCR shadow peak model includes 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 can be categorized into three cases: 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 follows:

[0089]

[0090] After obtaining the PCR shadow peak model, this method requires estimating the model parameters using sequencing reads of a given sample. For a given individual i, n are located at STR loci. i A read that completely spans the locus will cause the PCR peaks to aggregate as follows: in Let represent the difference between the observed allele and the true allele in the j-th READ traversing the current locus for individual i. Then, the likelihood function for n individuals in the second-generation sequencing sample dataset 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 from N samples that completely span the current STR locus;

[0096] m1 indicates the number of repeat units in the reads that are one more than the actual alleles;

[0097] n2 represents the number of alleles in the reads that have one less repeat unit compared to the actual alleles;

[0098] S5-3. The maximum likelihood estimates of parameters u and d can be obtained through formulas (3) and (4):

[0099]

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

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

[0102] The allele frequencies of STR loci have a relatively stable distribution in the population, which plays an important role in estimating the posterior probability during STR genotyping. Especially in low-coverage sequencing data, accurate allele frequency distribution can greatly improve the accuracy of STR genotyping results. For some common STR loci, such as forensic STR loci, the allele frequencies in the population 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 readily available allele frequency distribution data online. Therefore, this invention can automatically estimate the allele frequency distribution of user-input sequencing sample data, and can continuously correct the locus allele frequency data as the sample size increases, thereby improving the accuracy of STR allele genotyping. In order to obtain accurate and reliable STR locus allele frequency distribution, this 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) candidate alleles must appear in at least 2 samples. For datasets where the sequencing samples are insufficient to estimate the allele frequencies of custom STR loci, this invention assigns the same frequency value to each potential STR allele, thereby eliminating the influence of prior probability on posterior probability estimation.

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

[0104] Because only a small subset of reads in low-coverage sequencing data can completely traverse the entire STR locus, and considering the influence of PCR peaks, it is difficult to accurately infer the alleles of a sample solely from the number of supporting reads. However, when the PCR peak ratio of the STR locus and the frequency of alleles in the population are known, this prior knowledge can be used to estimate the posterior probability of alleles, thereby significantly improving the accuracy of allele genotyping. This invention uses the estimated PCR 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 highest posterior probability is taken as the final STR genotyping result. The maximum posterior probability estimation process is as follows:

[0105] For a given sample, assuming it has n reads that completely span the current STR locus, the set of alleles formed by these reads can be represented as A = {a1, a2, ..., a...} nThen, the maximum posterior probability for each possible candidate allele a (a∈A) can be calculated using the following formula:

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

[0107]

[0108] f in formula (6) a Frequency of candidate allele a at STR locus 1; a i This represents the i-th observed allele at the current locus in the sample; a+1 and a-1 represent shadow alleles that have one more repeat unit and one less repeat unit compared to the given candidate allele a, respectively; others represent the observed allele a. i The number of repeating units is greater or less than that of a given candidate allele 'a'. Since autosomal STR loci have two alleles, this method assumes that the two alleles at the STR locus are j (j∈A) and k (k∈A) to estimate the posterior probability of this type of allele; thus, the joint probability of all pairwise allele combinations can be estimated:

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

[0110]

[0111] Test example:

[0112] 1. This invention was benchmarked against three mainstream STR genotyping software programs (HipSTR, Strait Razor, and STRINGS). The test data consisted of 3000 simulated WGS (single-end, 150-bp) sequencing data points with sequencing depths of 10× (1000), 30× (1000), and 50× (1000), respectively. The tested 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, HalimuretiSimayijiang, and Claus "STRinNGS v2.0: Improved tool for analysis and reporting of STR sequencing data." Forensic Science International: Genetics 48(2020):102331."

[0116] Test results show that the genotyping method of this invention achieves a 100% detection rate and a genotyping accuracy of 99.37% 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 low-coverage WGS sequencing data, this invention can also be applied to high-coverage targeted sequencing data or amplicon sequencing data. To verify its genotyping effect, this invention used 687 real targeted sequencing data (single-end, 400-bp), with an average sequencing depth of 1570×, and tested 53 commonly used STR loci.

[0118] Test results show that the detection rate and accuracy of the typing method of this invention are as high as 99.64% and 99.9%, respectively. In comparison, the detection rate 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. This invention greatly improves the STR typing speed through algorithm optimization. Using 10 high-coverage real sequencing data (average 3216×) for typing speed testing, this invention can complete STR allele typing in only 136 seconds, which is 79 times faster than HipSTR and more than 10,000 times faster than STRinNGS.

[0120] In summary, the STR genotyping method of this invention, specifically applied to second-generation sequencing data, has higher detection and genotyping 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 genotyping 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, and is not intended to limit the scope of protection of the present invention. Simple modifications or equivalent substitutions made by those skilled in the art to the technical solution of the present invention do not depart from the essence and scope of the technical solution of the present invention.

Claims

1. A method for STR typing of next-generation sequencing data, characterized in that, Includes the following steps: S1: Custom STR locus configuration file; S2: Constructing a k-mer index for flanking sequences of STR loci: Constructing 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 described in step S2, execute steps S4-S7 respectively; S4: Extract STR locus allele sequences: Extract STR allele sequences from each READ of the sequence alignment file by combining CIGAR value and k-mer index, and generate the allele set of STR loci; S5: Construct a PCR shadow peak model for STR loci, and estimate the parameters of the PCR shadow peak model using a second-generation sequencing sample dataset and the maximum likelihood estimation method. The PCR shadow peak model includes parameters. and , This represents the probability of having one more repeating unit compared to the true allele, while This indicates the probability of having one less repeat unit compared to the true allele; the distribution of the PCR shadow peak model is shown below: For a given STR locus, This indicates the number of repeating units that differ between the observed allele and the true allele, and is divided into three cases: Identical to true alleles ; It has one more repeat unit than the true allele. ; And one less repeating unit compared to the true allele. ; Therefore, the distribution of the PCR shadow peak model can be written out: ； S6: Estimating candidate allele frequencies at STR loci: Estimating candidate allele frequencies using second-generation sequencing sample datasets; S7: Infer the most likely alleles and their posterior probabilities for STR loci in the sequence alignment file: Calculate the posterior probabilities of all candidate alleles using the allele set extracted in S4, and infer the most likely alleles accordingly.

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 a 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 construction of the STR locus flanking sequence k-mer index adopts a dynamic k-mer index construction strategy, that is, selecting the optimal length based on the repetitiveness of the STR locus flanking sequences. The length of the k-mer subsequence at this locus is then used as the k-mer subsequence length. The relative position of k-mer in the flanking sequence Stored in a hash table to prepare for subsequent allele extraction.

5. The STR typing method according to claim 1, characterized in that, In S4, the specific operation for extracting the STR allele sequence is as follows: S4-1. Using the sequence alignment file as input, first calculate the flanking sequence position of the STR using the CIGAR value of each READ and the position of the corresponding STR site on the reference genome; S4-2. Then, the calculated READ flanking sequence is verified using the flanking sequence in the configuration file, and the STR allele is directly calculated in this way. S4-3. When an insertion or deletion in the CIGAR value occurs at the junction of the flanking sequence and the core sequence and the number of bases does not exceed twice the length of the STR repeat unit, the insertion or deletion in the 5' and 3' flanking sequences that are identical to the repeat unit and have an alignment error is moved back to the STR core region using the CIGAR value correction method, and finally the STR allele sequence is extracted using the corrected CIGAR value meter. S4-4. When the allele length of a sample on the STR genome exceeds 40% of the allele length of the reference genome at that locus, sequence alignment software often introduces completely erroneous CIGAR values. To compensate for the defects of the CIGAR method, the k-mer index method is used to extract alleles: for each READ sequence, the corresponding STR locus 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 specific steps for estimating the parameters of the PCR shadow peak model using the second-generation sequencing sample dataset and the maximum likelihood estimation method are as follows: S5-1, For a given individual At the current STR locus, A read that completely spans the locus will cause the PCR peaks to aggregate as follows: ,in This represents the difference between the observed allele and the true allele in the j-th READ that spans the current locus for individual i. Then, for the second-generation sequencing sample dataset The likelihood functions for each individual are shown below: ； S5-2. Based on formulas (1) and (2), we can calculate the following about and Likelihood function: ； ； in: express The sum of the number of reads that completely span the current STR locus in each sample; This indicates the number of repeat units in the reads that are one more than the actual alleles; This indicates that the number of repeat units in the reads is one less than the number of true alleles; S5-3. The parameters can be estimated using formulas (3) and (4). and Maximum likelihood estimate: 。 7. The STR typing method according to claim 1, characterized in that, In S6, the estimation of candidate allele frequencies specifically involves: statistically analyzing the allele frequency distribution of the current STR locus using the second-generation sequencing sample dataset input by the user; and re-statistically analyzing the candidate allele frequencies of the current locus when the sample size increases.

8. The STR typing method according to claim 1, characterized in that, In S7, the most probable allele and its posterior probability for the current STR locus in the sequence alignment file are inferred. The specific operation is as follows: S7-1. If the current STR locus is defined as haplotype in the configuration file, it has [a certain number of STR loci]. If there are reads that completely span the locus, then the set of alleles observed in these reads can be represented as follows: Then, for each possible candidate allele The posterior probability can be calculated using the following formula: ； in, Indicates candidate alleles Frequency in the current STR locus; This indicates the sample's position at the current locus. One observed allele; and They represent the corresponding candidate alleles. Compared to the shadow peak alleles with an added repeat unit and those with a removed repeat unit; "Others" represents the observed alleles. Compared to a given candidate allele The number of repeating units may be greater or less than a certain number; 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, the two alleles of the STR locus are respectively... and k,j , k Therefore, the joint posterior probability of all pairwise allele combinations can be estimated: ； in, Indicates candidate alleles and alleles Frequency in the current STR locus; This indicates the sample's position at the current locus. One observed allele; and They represent the corresponding candidate alleles. Compared to the shadow peak alleles with an added repeat unit and those with a removed repeat unit; S7-3. By calculation, the posterior probabilities of all candidate alleles or allele combinations at the current STR locus are obtained. The allele or allele combination with the highest posterior probability is selected as the most likely allele or allele combination in the sequence alignment file at the current STR locus.

Images