STR automatic typing and naming analysis method and system based on high-throughput sequencing technology, terminal equipment and storage medium

Abstract

The invention belongs to the technical field of biological information, and particularly relates to an STR automatic typing and naming analysis method and system based on a high-throughput sequencing technology, terminal equipment and a storage medium. The analysis method comprises the following steps: setting an analysis configuration file through an STR flanking sequence and an STR repeated structure form of an STR locus to be detected, obtaining sequence text information of an hg38 reference genome and sequence text information of sequencing data according to the analysis configuration file, and comparing the sequence text information of the hg38 reference genome and the sequence text information of the sequencing data to obtain a core region of each STR locus of the sequencing data; therefore, the STR gene loci are automatically typed and named. A data analysis method which is not limited to a sequencing platform, simple to operate and visual in analysis result is established, the STR gene loci can be subjected to typing based on length polymorphism and sequence polymorphism, naming conforming to the international forensic genetics (ISFG) standard is performed, and flanking sequences of the STR gene loci can be analyzed at the same time.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics technology, and particularly relates to an analysis method, a system, a terminal device and a storage medium for automatic STR genotyping and naming based on high-throughput sequencing technology. Background Art

[0002] The short tandem repeat (STR) genotyping technology is a commonly used detection method in forensic physical evidence. Currently, the most widely used method is based on the polymerase chain reaction-capillary electrophoresis (PCR-CE) technology platform. This method only detects the number of repeats in the core region of the STR locus, and those with the same number of repeats are counted as the same allele. In recent years, with the development of high-throughput sequencing technology, the detection of STR loci has expanded from only focusing on length polymorphism to sequence polymorphism.

[0003] Currently, relevant commercial kits have been developed, such as the Forenseq TM DNA Signature Prep Kit and PowerSeq TM systems for the MiSeqFGx Forensic Genomics sequencing platform, and the Precision ID Globalfiler TM system for the IonTorrent platform, etc. For the analysis and reading of these sequencing results, corresponding data analysis tools have been developed, such as Universal Analysis Software (UAS), Torrent Suite Software (TSS), and Converge, etc. These tools have the advantages of simple operation and user-friendliness.

[0004] However, due to the closed software environment, users can only analyze the preset STR locus information, and the analysis of the flanking sequences of STR loci is limited. Users cannot modify various parameters to achieve personalized applications, which greatly hinders the development of high-throughput sequencing of STR loci and the accumulation of population data. At the same time, these software all need to be paid for use, which also poses an obstacle to the popularization of high-throughput sequencing technology in grass-roots forensic physical evidence laboratories. In addition, the genotyping of high-throughput sequencing STR requires a high level of bioinformatics knowledge reserve for users, so it has also become a "stumbling block" problem for most grass-roots units and researchers. Summary of the Invention

[0005] In view of the above deficiencies, the present invention provides an analysis method for automatic genotyping and naming of STR based on high-throughput sequencing technology. The present invention establishes a data analysis method that is not limited to the sequencing platform, has simple operations, and intuitive analysis results. It can genotype STR loci based on length polymorphism and sequence polymorphism, perform naming in line with the standards of the International Society for Forensic Genetics (ISFG), and simultaneously analyze the flanking sequences of STR loci.

[0006] The technical solution of the present invention includes:

[0007] In a first aspect, the present invention provides an analysis method for automatic genotyping and naming of STR based on high-throughput sequencing technology, characterized in that the analysis method includes the following steps:

[0008] S1. Determine the flanking sequences and repeat structure forms of the core regions of the STR loci to be tested based on the hg38 reference genome, referred to as STR flanking sequences and STR repeat structure forms;

[0009] S2. Set an analysis configuration file according to the STR flanking sequences and STR repeat structure forms;

[0010] S3. Establish an hg38 index for the hg38 reference genome, and obtain the text information of the STR repeat structure sequences and flanking sequences of the hg38 reference genome according to the analysis configuration file in step S2, referred to as reference genome text information;

[0011] S4. After aligning the sequencing data with the hg38 reference genome, obtain the sequencing data sequence text information according to the analysis configuration file in step S2;

[0012] S5. Compare the reference genome text information in step S3 with the sequencing data sequence text information in step S4, so that the flanking sequences of the hg38 reference genome are anchored to the STR flanking sequences, and the core sequences of each STR locus in the sequencing data are located in the middle position, referred to as the core regions of each STR locus in the sequencing data;

[0013] S6. Name the core regions of each STR locus in the sequencing data in step S5 according to the repeat structure form obtained in step S1;

[0014] S7. Run the analysis command to obtain the sequenced data file; the analysis command is: pythonSTRanalysis.py -i example.fastq.gz -o / A -n example.ini -r hg38.fa -t 50 -f; where, -i represents the input high-throughput sequenced data file to be analyzed; -o represents the location to define the output file; -n represents the analysis configuration file described in step S2; -r represents the hg38 reference genome; -t represents the number of threads specified for STR analysis; -f represents the forced execution of the analysis.

[0015] The sequenced data file includes: the naming and repeat structure forms of the core regions of each STR locus in the sequencing data.

[0016] S8. Summarize the reads information belonging to the same STR locus in the sequenced data file. There are various repeat structure forms in the core regions of each STR locus in the sequencing data. Count each repeat structure form and determine the genotype based on the count of the repeat structure form.

[0017] Specifically, the method for determining the core region of the STR locus to be tested described in step S1 is: in the hg38 reference genome, find the base sequence position containing the STR locus, which is the core region of the STR locus to be tested.

[0018] Specifically, the method for determining the STR flanking sequences described in step S1 is: the upstream sequence and downstream sequence of the core region of the STR locus to be tested are the upstream flanking sequence and downstream flanking sequence respectively, which are the STR flanking sequences.

[0019] Preferably, the length of the STR flanking sequences is 20 - 30 bases.

[0020] Specifically, the method for determining the STR repeat structure form described in step S1 is: find the base repeat sequence in the core region of the STR locus to be tested; calculate the number of repeats.

[0021] Preferably, the base repeat sequence is a base fragment with a length of 3 - 6 bases that continuously repeats 2 times or more in the core region of the STR locus to be tested;

[0022] Specifically, the manifestation forms of the STR repeat structure form described in step S1 include:

[0023] a. When the STR locus contains only one base repeat sequence, the manifestation form of the STR repeat structure form is: STR locus [X1]n;

[0024] Where, X1 represents the base repeat sequence; n represents the number of repeats, which is a positive integer greater than 1.

[0025] b. When the STR locus contains y types of base repeat sequences, the manifestation form of the STR repeat structure is: STR locus [X1]n[X2]n…[Xy-1]n[Xy]n;

[0026] Among them, y represents the number of types of repeat sequences, which is a positive integer greater than 1; X1 represents the first type of base repeat sequence from the 5'-end to the 3'-end of the core region of the STR locus to be detected; X2 represents the second type of base repeat sequence from the 5'-end to the 3'-end of the core region of the STR locus to be detected; Xy-1 represents the second-to-last type of base repeat sequence from the 5'-end to the 3'-end of the core region of the STR locus to be detected; Xy represents the last type of base repeat sequence from the 5'-end to the 3'-end of the core region of the STR locus to be detected; n represents the number of repeats, which is selected from positive integers greater than 1 that are the same or different;

[0027] c. When there are non-repeat base fragments between the same base repeat sequences or different base repeat sequences, the manifestation forms of the non-repeat base fragments include: Z[1] or Nz; among them, Z represents a single base; z represents the number of bases, which is selected from positive integers greater than 1.

[0028] When the length of the non-repeat base fragment is 1 base, the manifestation form of the non-repeat base fragment is Z[1]; the Z is selected from: adenine, thymine, guanine or cytosine; when the length of the non-repeat base fragment is more than 1 base, the manifestation form of the non-repeat base fragment is Nz.

[0029] Specifically, the analysis configuration file in step S2 includes necessary parameters; the necessary parameters include: STR locus name, chromosome where the core region of the STR locus is located, chromosomal position where the core region sequence of the STR locus starts, chromosomal position where the core region sequence of the STR locus ends, upstream flanking sequence length, downstream flanking sequence length, number of tolerated bases in the upstream flanking sequence, number of tolerated bases in the downstream flanking sequence, number of bases of the base repeat sequence in the core region of the STR locus, naming method of the core region of the STR locus; the naming method of the core region of the STR locus is confirmed according to the STR repeat structure form.

[0030] Specifically, the tool used to establish the hg38 index in step S3 is the samtools tool; the format of the hg38 index is the fai format; the implementation command to establish the hg38 index is: samtools faidx hg38.fa.

[0031] Specifically, a flanking sequence tolerance mechanism is established during the alignment process in step S5.

[0032] Preferably, the flanking sequence error tolerance mechanism includes:

[0033] If a completely matched flanking information is identified in the text information of the sequencing data sequence, the core sequence of the STR locus obtained by reads is used; if a completely matched flanking information is not identified in the text information of the sequencing data sequence, the number of error tolerance bases of the flanking sequence is increased, and the reads with no more than 2 mismatched bases are classified into the STR locus.

[0034] Specifically, the determination method for counting and judging the genotype according to the repeated structure form in step S8 is as follows:

[0035] If the repeated structure form with the most counts exceeds 70% of the reads of all the repeated structure forms of this locus, it is judged that this locus is homozygous; if the repeated structure form with the most counts does not exceed 70% of all the reads, the sum of the reads of the most and the second most repeated structure forms is calculated. If it exceeds 70% of all the reads, it is judged that this locus is heterozygous.

[0036] In a second aspect, the present invention provides an STR automatic genotyping and naming analysis system based on high-throughput sequencing technology, and the analysis system implements the above analysis method.

[0037] In a third aspect, the present invention provides a terminal device, which includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the above analysis method is implemented.

[0038] In a fourth aspect, the present invention provides a computer-readable storage medium, which stores a computer program. When the computer program is executed by a processor, the above analysis method is implemented.

[0039] The beneficial effects of the present invention are as follows:

[0040] (1) The analysis method of the present invention can analyze any STR locus of interest and can cope with various applications by changing various parameters, getting rid of the limitation that commercial software can only analyze inherent STRs. The application scope of this method is not limited to forensic applications. This method can be compatible with sequencing data of poor quality and has the ability to detect rare mutations.

[0041] (2) The analysis method of the present invention can directly output quantitative information such as heterozygosity and the balance of multiplex amplification sites in the results, and users can directly judge the credibility of the results of this locus according to the values.

[0042] (3) In addition to forensic applications, the analysis method of the present invention can also be used for STR analysis for various other purposes. By establishing a fault tolerance mechanism for fuzzy alignment of flanking sequences, sufficient information mining is ensured. At the same time, functions such as automatic naming of STR loci based on sequence polymorphism, calculation of population genetics data, automatic judgment of homozygosity / heterozygosity, automatic identification of stutter artifacts, and detection and prompt of rare mutations can be completed. The operation is simple and the running speed is fast. It has been verified that more than 150 loci can be detected simultaneously, and the detection results do not affect each other.

[0043] (4) The analysis method of the present invention can enhance the probative force contained in a single STR locus because the results obtained by this analysis method simultaneously reflect length polymorphism and sequence polymorphism, providing more polymorphism information than traditional methods that only show length polymorphism. In the analysis that has been carried out, based on the analysis of the length polymorphism of 65 autosomes, the cumulative probability of individual identification is 1 - 4.62153E - 70, while based on sequence polymorphism, the cumulative probability of individual identification is 1 - 8.98372E - 76, with a difference of approximately 510,000 times. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 It is a schematic diagram of the STR locus D16S539.

[0045] Figure 2 It is the output file of the BGI sequencer G99.

[0046] Figure 3 It is the output file of the Illumina platform sequencer MiSeq FGx System.

[0047] Figure 4 It is the output file of the Element Biosciences sequencer element AVITI. DETAILED DESCRIPTION OF THE INVENTION

[0048] The present invention will be further described in detail below in conjunction with specific embodiments. The following embodiments are not used to limit the present invention, but only to illustrate the present invention. The experimental methods used in the following embodiments are not specifically described, and the experimental methods without specific conditions in the embodiments usually follow conventional conditions. The materials, reagents, etc. used in the following embodiments can be obtained from commercial channels without specific description.

[0049] Example 1 Analysis Environment Setup

[0050] Install the Python running environment. The normal operation of the analysis method of the present invention requires a support environment and software. These environments and software have been packaged and can be installed and used at any time. Once installed, it can be used for a long time. Versions above Python 3.x are currently widely used Python tools. At the same time, in the analysis environment, it is necessary to install the samtools software for building indexes, the hisats2 tool for text information comparison, and the agrep tool for fuzzy retrieval of text. This analysis environment can be installed and downloaded in a packaged manner.

[0051] Example 2 Determination of the flanking sequence and repeated structure form of the core region and setting of the analysis configuration file

[0052] 1. Method for determining the flanking sequence and repeated structure form of the core region of STR loci

[0053] According to the STR loci to be detected, on the hg38 reference genome, determine the upstream and downstream sequences of the core repeat region of the STR to be detected, that is, the "flanking sequences". To ensure the uniqueness of the detection of the STR core repeat region, the upstream and downstream flanking sequences are approximately 25 bp; at the same time, based on the hg38 reference genome, set the possible repeated structure forms of the STR for subsequent automatic naming.

[0054] Problem of the selection of flanking length: Determining the flanking sequence is a necessary condition for STR analysis. If the flanking sequence is selected too short, combined with the error tolerance mechanism of the analysis method of the present invention, it is more likely to match multiple incorrect positions on the reference genome; at the same time, if the flanking sequence is selected too long, first, the analysis speed will be slowed down, and at the same time, high requirements are also put forward for the sequencing quality, which is not suitable for the analysis of samples with poor sequencing quality. Therefore, according to experience, the flanking sequences used in the analysis method of the present invention are approximately 20 - 30 bases.

[0055] Examples of confirming the flanking sequence and repeated structure form of the core region of STR loci:

[0056] (1) Taking the D10S1248 locus as an example of the STR locus to be detected: In the hg38 reference genome, the base sequence containing this locus is chr10:129294215 - 129294323, and the sequence is SEQ ID NO.1:

[0057] ATATTAATGAATTGAACAAATGAGTGAGT GGAAGGAAGGAAGGAAG GAAGGAAGGAAGGAAGGAAGG AAGGAAGGAAGGAA ATGAAGACAATAC AACCAGAGTTGTTCC.

[0058] In this segment of the sequence, a repeating structure of the four-base [GGAA] can be observed. The specific repeating positions are chr10:129294244 - 129294295 (the underlined part in SEQ ID NO.1), with a total of 13 repetitions. Then this position is defined as the STR core repeating sequence, simply referred to as the core region. The two ends of this sequence, namely chr10:129294215 - 129294243 and chr10:129294296 - 129294323, are the upstream and downstream flanking sequences respectively. The above-mentioned repeating position is called the D10S1248 locus, denoted as D10S1248[GGAA]13, which means that at this locus in the reference genome, the GGAA four-base structure repeats 13 times.

[0059] D10S1248 has good polymorphism in the population. Regardless of whether the repeating situations of the studied samples at this locus are the same, the flanking sequences use the same chromosomal positions (that is, the upstream flanking sequence ends at the position chr10:129294243, and the downstream flanking sequence starts at the position chr10:129294296), which can be simply described as "clamping the two ends of the flanks and calculating the number of core repetitions". The core sequence of this locus consists only of the four-base GGAA repeats, and the repeating structure form is set as: "D10S1248[GGAA]n".

[0060] The number of repeating bases in this core region is 4 bases, and this value is used for the automatic operation of the STR length-based naming in the analysis method. n represents any positive integer that may occur. During the analysis of the sample, the analysis method of the present invention can distinguish the number of bases at this locus of the sample. Dividing the number of bases by 4 gives the n value. If it is found during the analysis that the core sequence does not conform to the preset form (such as the length is not an integer multiple of 4, or does not strictly follow the GGAA repeating structure), a prompt will be given in the analysis result. This situation occurs occasionally, which is caused by the mutation of individuals, and this further reflects the sensitivity of the analysis method of the present invention.

[0061] (2) Taking the TH01 locus as an example of the STR locus to be measured: In the hg38 reference genome, the base sequence containing this locus is chr11:2171065 - 2171134, and the sequence is SEQ ID NO.2:

[0062] AGGGAACACAGACTCCATGGTG AATGAATGAATGAATGAATGAATG AATG AGGGAAATAAGGGAGGAAC.

[0063] In this segment of the sequence, a repeating structure of the four-base sequence [AATG] can be observed. The specific repeating positions are chr11:2171088 - 2171115 (the underlined part in SEQ ID NO.2), with a total of 7 repetitions. Then this position is defined as the STR core repeating sequence. The two ends of this sequence, namely chr11:2171065 - 2171087 and chr11:2171116 - 2171134, are the upstream and downstream flanking sequences respectively. This locus is denoted as TH01

[0064] [AATG]7, with a length of 28 bases. However, through the study of more than a hundred unrelated individuals, other forms were found in the core region of this locus. This situation is caused by mutations in the population. For example, the TH01 [AATG]5ATG AATG form, that is, [AATG] repeats 5 times, followed by ATG once, and then [AATG] once, with a length of 27 bases. Then the repeating structure form of TH01 is set as: "TH01[AATG]n" or "TH01[AATG]n ATG[1][AATG]n".

[0065] When the repeating structure form is: TH01[AATG]n ATG[1][AATG]n, it means that [AATG] appears any positive integer number of times, then ATG appears a definite 1 time, and then [AATG] appears any positive integer number of times. If one of the positive integers n is found to be 1 in the analysis, then in the result file, it will directly appear in the form of AATG without the "[]" symbol. If the analysis method of the present invention recognizes this form, it will adjust the naming form based on the base length. For example, the digital naming of TH01[AATG]5ATG AATG is TH01 6.3. The number before the decimal point represents the number of bases satisfying 6 (5 + 1) repeating units (i.e., 24), and the number after the decimal point represents the number of extra bases when satisfying 6 repeating units but less than 7 repeating units. This is a well-established naming form in the art.

[0066] (3) Taking the D12ATA63 locus as an example of the STR locus to be measured: In the hg38 reference genome, the base sequence containing this locus is chr12:107928571 - 107928653, and the sequence is SEQ ID NO.3:

[0067] GGATAGCAATTTAAAAATG TTGTTGTTGTTATTATTATTATTATTATTA TTATTATTA CTTGAGACAGGGTCTCGCTCTGTTA。

[0068] In this segment of the sequence, the repetition of the [TTG] and [TTA] sequences can be observed. The specific location is chr12:107928590-107928628 (the underlined part in SEQ ID NO.3), which contains 39 bases and is presented as D12ATA63[TTG]3[TTA]10. This locus is denoted as D12ATA63[TTG]3

[0069] [TTA]10, and the repeat structure form is set as: "D12ATA63[TTG]n[TTA]n". This locus contains two or more core sequences and is called a complex STR. When performing length-based naming, the entire STR core region should be calculated together.

[0070] (4) Taking the D14S608 locus as an example of the STR locus to be tested: In the hg38 reference genome, the base sequence containing this locus is chr14:28380240-28380359, and the sequence is SEQ ID NO.4:

[0071] GTGGTACAGGTAGATAAATGGAT GATAGATAGATA ATAGAGATAGAT GATAGACAGATAGATAGATA GATAGATAGATAGATAGATAGATA GAGTA TATATATAATGGACTATATAATAT.

[0072] In this segment of the sequence, the repetition of the [GATA] sequence can be observed. The specific locations are chr14:28380263-28380274 and chr14:28380287-28380330 (the underlined part in SEQ ID NO.4). Between the two core regions, there is an irregular sequence with a length of 12 bases that does not form a repeat structure. In actual analysis, this segment of the sequence is denoted as N12, representing 12 arbitrary bases. The repeat structure form of this locus is set as: "D14S608[GATA]n N

[12] [GATA]nGACA[GATA]n". Where n is a positive integer, and other key points are the same as the above situation.

[0073] 2. Set the analysis configuration file for the STR locus

[0074] In this embodiment, taking the STR locus D1S1677 as an example, the locus information included in the profile is shown, including the chromosomal position (start) where the core region sequence begins, the chromosomal position (stop) where the core sequence ends, the lengths of the upstream and downstream sequences analyzed (flank_up_length, flank_down_length), the number of tolerated base mismatches (mism_up, mism_down), the number of base pairs in the core sequence repeat unit (unit_length), the nomenclature method (nomenclature), and other necessary parameters. The necessary parameters are as follows:

[0075] “[D1S1677]

[0076] chr = chr1

[0077] start = 163590026

[0078] stop = 163590085

[0079] flank_up_length = 25

[0080] flank_down_length = 25

[0081] mism_up = 2

[0082] mism_down = 2

[0083] unit_length = 4

[0084] nomenclature = TTCC[+]”.

[0085] It shows that in the hg38 reference genome, the chromosomal position of the core region of D1S1677 is: chr1:163590026 - 163590085. In this example, the lengths of the upstream and downstream flanking sequences are both 25, that is, chr1:163590001 - 163590025 is used as the upstream flank, and chr1:163590086 - 163590110 is used as the downstream flank sequence. mism_up and mism_down are used as the number of tolerated mismatches for the upstream and downstream flanks respectively. The number of tolerated mismatches is beneficial for capturing more information in the case of poor sequencing quality. unit_length represents the number of repeated bases in the core region, and nomenclature represents the structural form of the STR core region.

[0086] Other optional parameters such as background noise filtering threshold, minimum number of reads threshold, heterozygosity threshold, etc. can further put forward personalized requirements for STR analysis. Information such as the lowest limit of the detection signal and the number of detectable alleles can be set, and such parameters are set based on experience. Other optional parameters are as follows:

[0087] “subtract=0.0

[0088] strand=f

[0089] noise_filter=0.01

[0090] min_reads=50

[0091] max_num_unique=6

[0092] min_frac_genotype=0.7

[0093] min_frac_profile=0.004

[0094] hetero_balance=0.25,0.75

[0095] max_reads_unique_not_called=0.15

[0096] min_unex=20

[0097] nomenclature=TTCC[+]”。

[0098] The format of this configuration file is in ini format and is named example.ini in this embodiment.

[0099] Example 3: Establish an index, perform analysis in a Python environment, and output the analysis results

[0100] 1. Use samtools to build an index for the hg38 reference genome

[0101] In the analysis profile described in Example 2, the chromosomal positions of the loci and flanking sequences are given, but the specific sequence information is not directly provided. When performing analysis, it is necessary to retrieve the reference genome based on the provided chromosomal positions and capture the base information at the corresponding positions. This process requires building an index. The index built by samtools is in the.fai format, which is a commonly used index form in gene analysis and is well-known information in the industry. For the reference genome selected for alignment analysis, generally the fasta format file of the latest version of the hg38 reference genome, namely hg38.fa, is used. The bioinformatics software samtools is used to build an index for the hg38.fa reference genome, and the generated index is hg38.fa.fai. Specific implementation command: samtools faidx hg38.fa.

[0102] 2. Obtaining the sequence text information of the sequencing data

[0103] According to the above content, the specific base information of the flanks can be obtained through the chromosomal positions of the flanks and based on the hg38 reference genome index. The flank position information is converted into base information, and this information is captured. The hisats2 tool is used to align the original demultiplexed data with the hg38 reference genome, and the sequence information is converted into text information. At this time, a large amount of repeated or non-repeated text information will be obtained, which is the sequence text information.

[0104] 3. Obtaining the STR repeat core sequence of the sequencing data

[0105] The T.R.E. mode of the agrep tool is used to perform a fuzzy alignment between the flank base information of the hg38 reference genome and the sequence text information in the sequencing data. The upstream and downstream of the relevant loci in the sequencing data are anchored to the upstream and downstream in the hg38 reference genome respectively, and the middle position is the STR repeat core sequence.

[0106] If a completely matching flank information can be identified, then this read (reads refers to the information in the demultiplexed file, which is commonly referred to as reads in the industry) can be matched to the corresponding STR locus. Then, the STR core region is further counted and named. It should be noted that to ensure the uniqueness of the match, the upstream and downstream flank sequences should ensure a unique match in the reference genome, and at the same time, in the demultiplexed file reads, the upstream and downstream flanks must be paired and matched to be counted in a certain STR locus.

[0107] During this process, a flanking sequence error tolerance mechanism can be established for fuzzy alignment, which can maximize the capture and summary of core sequence information. If the sequencing quality is poor and the flanking sequences in the reads cannot be fully matched with the flanking sequences determined by the reference genome, the error tolerance mechanism will come into play at this time, and the reads with no more than 2 base mismatches can be classified into this STR locus (when the number of mismatches is set to 2).

[0108] Explanation of the fuzzy alignment error tolerance mechanism: The error tolerance mechanism can maximize the role of data mining when the sequencing quality is poor. However, when the number of error tolerance bases increases, the analysis speed will slow down. If the specific sequences detected by the error tolerance mechanism exceed the threshold, the reads containing this specific mutation will be summarized and output, prompting the user of the detection of rare mutations. The following example demonstrates the working principle of the error tolerance mechanism:

[0109] As Figure 1 shown, this figure shows a schematic diagram of the STR locus D16S539. In the hg38 reference genome, the number of repeats of this STR is 11 times. The following are examples of sequencing reads respectively. During detection, when the flanking error tolerance base number is 0, a total of 1000 reads can be detected. When the error tolerance is 1, a total of 5000 reads can be detected, and so on. When the error tolerance is 2, 10000 reads can be detected, and when the error tolerance is 3, 10500 reads can be detected. The upstream and downstream flanking error tolerance numbers can be set independently without affecting each other, maximizing the detection ability of low-quality sequencing data.

[0110] The same sample is analyzed using different flanking error tolerance numbers, and the results are shown in Table 1:

[0111] Table 1 Analysis results of the same sample using different flanking error tolerance numbers

[0112]

[0113] When the flanking error tolerance number is smaller, the analysis time-consuming is shorter, and the simultaneous analysis of 5 loci can be completed within a few minutes. This analysis method only needs to input one analysis code, with simple operation, intuitive results, and fast analysis speed.

[0114] 4. STR Typing and Naming

[0115] For the obtained STR repeat core sequence, according to the number of bases in the STR repeat unit, the length information is automatically calculated, and according to the possible repeat structural forms of the STR obtained in Example 2, the locus is automatically named. At the same time, the obtained STR repeat sequences are counted, and according to the default parameters or set parameters, it is automatically analyzed whether this locus is homozygous or heterozygous, the stutter pseudo-signal is automatically discriminated, and at the same time, according to the count, the heterozygosity and balance information are automatically output.

[0116] Run the analysis command in the configured Python environment (i.e., the automatically installed and configured environment in Example 1): python STRanalysis.py -i example.fastq.gz -o / A -n example.ini -r hg38.fa -t 50 -f. The character descriptions are as follows:

[0117] -i: Input the high-throughput sequencing demultiplexed file to be analyzed; -o: Define the location of the output file; -n: Specify the.ini format analysis configuration file used for STR analysis; -r: Specify the reference genome sequence for STR analysis. An.fai index file must be established for the reference genome; -t: Specify the number of threads for STR analysis. The more threads, the faster the analysis speed; -f: Set the forced execution of the analysis. Using this parameter can overwrite the previous analysis results and only retain the new analysis results, but the log files of previous analyses will be retained.

[0118] 5. Detection of flanking special sequences (optional)

[0119] The analysis method of the present invention can perform special processing on flanking positions of interest (such as SNP sites) or avoidance positions (such as indel sites). For positions of interest, the base sequence of this position can be separately implemented without occupying the flanking tolerance base number; for avoidance positions, they can be ignored during analysis to avoid misreading of STR caused by flanking sequence shift.

[0120] 6. Output analysis results

[0121] Input the command in step "2.3 STR genotyping and naming" of this example and execute it. The analysis system will check whether the code statements are correct and whether the called files are available. If the check is correct, it will perform analysis and output a csv format file for direct viewing; if an error is detected in one of the files, an error message will be prompted and the analysis will be aborted for further modification.

[0122] Summarize the reads information of the same STR locus in the sequencing demultiplexed file. There may be multiple repeat forms in the core region, and count according to each repeat form. If the one with the most counts exceeds 70% of the total reads of this locus, it is judged that this locus is homozygous; if the one with the most counts does not exceed 70% of the total reads, calculate the sum of the reads with the most and the second most counts. If it exceeds 70% of the total reads, it is judged that this locus is heterozygous.

[0123] Example 3 The method for STR automatic genotyping and naming is used to test the demultiplexed files of multiple sequencing platforms

[0124] The analysis method of the present invention is applicable to the off-machine files of multiple sequencing platforms, and has been tested on the Illumina sequencing platform, the BGI sequencing platform, and the Element AVITI sequencing platform, and the results are reliable.

[0125] Test 1: STR analysis of the off-machine file of the BGI sequencer G99. The following is a screenshot of the off-machine file in fastq.gz format. From the identifier information line (the line starting with @), it can be seen that this off-machine file comes from the BGI sequencing platform. This off-machine file is single-end 400-read sequencing (SE 400) on the BGI sequencing platform, see Figure 2 .

[0126] Simultaneously analyze the D1S1656, D2S1338, D3S1358, vWA, and Penta D loci. Analysis running code: "python STRanalysis.py -i FT100012298_L01_17-25.fq.gz -o huadaanalysis-nexample.ini -r hg38.fa -t 50". The analysis result file results.csv is displayed in the output folder huadaanalysis, and the information can be directly viewed. The analysis results are shown in Table 2:

[0127] Table 2 Analysis results of Test 1

[0128]

[0129]

[0130]

[0131]

[0132]

[0133] The above table shows the name of the STR locus (Locus), the core sequence length (Length), the detected alleles (Allele), the alleles determined to be real (PredAllele), whether it is the same as the repeat structure form (ExactParse), the reads count (count), the signal ratio of the allele at this locus (read_fraction), the description of unexpected results (Comment), the STR core sequence bases (STR), the upstream flanking sequence bases (flank_up_match) and the number of tolerated differences (up_fl_num_diffs), the downstream flanking sequence bases (flank_down_match) and the number of tolerated differences (down_fl_num_diffs), the real composition form of the core sequence (ReportedNom), the heterozygosity information (Hb), the linkage equilibrium information (Lb), and the STR locus naming based on sequence information (STRidER) given according to the naming specification. As can be seen from Table 2, this sample is heterozygous at all 5 loci in example.ini.

[0134] Test 2: STR analysis of the demultiplexed file of the Illumina platform sequencer MiSeqFGx System. The following is a screenshot of the demultiplexed file in fastq.gz format. From the identifier information line (the line starting with @), it can be seen that this demultiplexed file is from the Illumina sequencing platform. This demultiplexed file was sequenced with paired-end 300 reads (PE300) on the Illumina sequencing platform. See Figure 3 。

[0135] Analyze the loci in example.ini simultaneously and run the following code: "python STRanlysis.py -i 100_S10_L001_R1_001.fastq.gz -o illuminaanalysis -n example.ini -r hg38.fa -t 50". The analysis result file results.csv is displayed in the output folder illuminaanalysis, and the information can be directly viewed. The analysis results are shown in Table 3:

[0136] Table 3 Analysis Results of Test 2

[0137]

[0138]

[0139]

[0140]

[0141] As can be seen from the information in Table 3, the sample is homozygous at the D1S1656 locus and heterozygous at the remaining loci.

[0142] Test 3: STR analysis of the output file of the Element Biosciences sequencer element AVITI. The following is a screenshot of the output file in fastq.gz format. From the identifier information line (the line starting with @), it can be seen that this output file is from the element sequencing platform. This output file was sequenced with paired-end 300 reads (PE300) on the element sequencing platform. See Figure 4 .

[0143] Analyze the loci in example.ini simultaneously and run the following code: "python. / src / STRinNGS2.py -i 92_R1.fastq.gz -o elementanalysis -n example.ini -r hg38.fa -t 50". The analysis result file results.csv is displayed in the output folder elementanalysis, and the information can be directly viewed. The analysis results are shown in Table 4:

[0144] Table 4 Analysis Results of Test 3

[0145]

[0146]

[0147]

[0148]

[0149]

[0150] The analysis results in Table 4 show rich STR information. If there is too much redundant or complex information in the sequencing results, the results.csv file also gives a description and judgment of the situation (in the comment column).

[0151] To demonstrate the universality of this analysis method across sequencing platforms, when analyzing the above three samples, exactly the same example.ini configuration file was used, including single-end and paired-end sequencing, different sequencing read lengths, and different sequencing platforms, which proves the convenience and reliability of this analysis method. The detailed information of each sequence is shown in Table 5:

[0152] Table 5 Detailed Information of Each Sequence

[0153]

[0154]

[0155]

[0156]

[0157]

[0158]

[0159]

[0160] The above detailed description is a specific description of one of the feasible embodiments of the present invention, and this embodiment is not intended to limit the patent scope of the present invention. The protection scope of the present invention shall be subject to the appended claims.

Claims

1. A method for automatic STR typing and naming based on high-throughput sequencing technology, characterized in that: The analytical method comprises the following steps: S1. Determine the flanking sequence and repeat structure form of the core region of the STR locus to be tested based on the hg38 reference genome, which are called STR flanking sequence and STR repeat structure form; S2, setting analysis profiles according to STR flanking sequences and STR repeat structure forms; S3, establishing an hg38 index for the hg38 reference genome, and obtaining text information of the STR repeat structure sequence and flanking sequence of the hg38 reference genome according to the analysis configuration file of step S2, which is referred to as reference genome text information; S4, after aligning the sequencing data with the hg38 reference genome, obtaining the sequencing data sequence text information according to the analysis configuration file of step S2; S5, comparing the reference genome text information of step S3 with the sequencing data sequence text information of step S4, anchoring the flanking sequence of the hg38 reference genome with the STR flanking sequence, and the middle position is the core sequence of each STR locus in the sequencing data, which is called the core region of each STR locus in the sequencing data; S6, naming the core region of each STR locus of the sequencing data of step S5 according to the repeating structure form obtained in step S1; S7, run the analysis command to obtain the sequencing file; the analysis command is: python STRanalysis.py -iexample.fastq.gz -o / An example.ini -r hg38.fa -t 50 -f; wherein -i represents the input of the high-throughput sequencing file to be analyzed; -o represents the location of the output file; -n represents the analysis configuration file described in step S2; -r represents the hg38 reference genome; -t represents the number of threads specified for STR analysis; -f represents the mandatory execution of the analysis; The sequencing download file includes: the naming and repeat structure of the core region of each STR locus in the sequencing data; S8. Summarize the reads information belonging to the same STR locus in the sequencing file. There are multiple repetitive structural forms in the core region of each STR locus in the sequencing data. Count each repetitive structural form and determine the typing based on the count of the repetitive structural form.

2. The analysis method according to claim 1, characterized in that The method for determining the core region of the STR locus to be tested in step S1 is as follows: in the hg38 reference genome, the base sequence position containing the STR locus is found, which is the core region of the STR locus to be tested; The method for determining the STR flanking sequence in step S1 is as follows: the upstream sequence and the downstream sequence of the core region of the STR locus to be tested are the upstream flanking sequence and the downstream flanking sequence, respectively, which are the STR flanking sequence; and the length of the STR flanking sequence is 20-30 bases.

3. The analysis method according to claim 1, characterized in that The method for determining the STR repeat structure form in step S1 is: finding the base repeat sequence in the core region of the STR locus to be tested; calculating the number of repetitions; the base repeat sequence is a base segment with a length of 3-6 bases that appears continuously for 2 or more times in the core region of the STR locus to be tested; The STR repeat structure described in step S1 may be expressed in the following forms: a. When the STR locus contains only one type of base repeat sequence, the STR repeat structure is expressed as: STR locus [X1]n; Wherein, X1 represents the base repeat sequence; n represents the number of repetitions, which is a positive integer greater than 1; b. When the STR locus contains y base repeat sequences, the STR repeat structure is expressed as: STR locus [X1]n[X2]n…[Xy-1]n[Xy]n; Wherein, y represents the number of repeating sequences, which is a positive integer > 1; X1 represents the first base repeating sequence from the 5' end to the 3' end of the core region of the STR locus to be tested; X2 represents the second base repeating sequence from the 5' end to the 3' end of the core region of the STR locus to be tested; Xy-1 represents the second to last base repeating sequence from the 5' end to the 3' end of the core region of the STR locus to be tested; Xy represents the last base repeating sequence from the 5' end to the 3' end of the core region of the STR locus to be tested; n represents the number of repetitions, which is selected from the same or different positive integers > 1; c. When there are non-repeating base segments between the same base repeating sequence or different base repeating sequences, the non-repeating base segments may be expressed as: Z[1] or Nz; wherein Z represents a single base; and z represents the number of bases, which is selected from a positive integer greater than 1; When the length of the non-repetitive base fragment is 1 base, the non-repetitive base fragment is expressed as Z[1]; the Z is selected from: adenine, thymine, guanine or cytosine; when the length of the non-repetitive base fragment is more than 1 base, the non-repetitive base fragment is expressed as Nz.

4. The analysis method according to claim 1, characterized in that The analysis configuration file described in step S2 includes necessary parameters; the necessary parameters include: STR locus name, chromosome where the STR locus core region is located, chromosome position where the STR locus core region sequence starts, chromosome position where the STR locus core region sequence ends, upstream flanking sequence length, downstream flanking sequence length, upstream flanking sequence error tolerance base number, downstream flanking sequence error tolerance base number, number of bases in the STR locus core region base repeat sequence, and STR locus core region naming method; the STR locus core region naming method is confirmed according to the STR repeat structure form.

5. The analysis method according to claim 1, characterized in that The tool used to establish the hg38 index in step S3 is the samtools tool; the format of the hg38 index is the fai format; the implementation command for establishing the hg38 index is: samtoolsfaidx hg38.fa.

6. The analysis method according to claim 1, characterized in that During the alignment process described in step S5, a flanking sequence error tolerance mechanism is established, and the flanking sequence error tolerance mechanism includes: If a fully matched flanking information is identified in the sequencing data sequence text information, the core sequence of the STR locus obtained by the reads is used; if a fully matched flanking information is not identified in the sequencing data sequence text information, the number of flanking sequence tolerance bases is increased, and the number of reads with less than 2 base mismatches is classified into the STR locus.

7. The method according to claim 1, characterized in that The method for determining the typing based on the count of the repeated structural forms described in step S8 is: If the repeat structure form with the most counts exceeds 70% of the total number of reads of the repeat structure form of the locus, the locus is judged to be homozygous; if the repeat structure form with the most counts does not exceed 70% of the total number of reads, the sum of the reads of the most and second most repeat structure forms is calculated. If it exceeds 70% of the total number of reads, the locus is judged to be heterozygous.

8. An analysis system for automatic STR typing and naming based on high-throughput sequencing technology, characterized in that: The analysis system implements the analysis method described in any one of claims 1-7.

9. A terminal device, characterized in that: The terminal device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, and the processor implements the analysis method described in any one of claims 1 to 7 when executing the computer program.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the analysis method according to any one of claims 1 to 7 is implemented.

Images