A method, electronic device and storage medium for detecting short tandem repeat expansion and genotyping

Abstract

This invention provides a method, electronic device, and storage medium for detecting short tandem repeat expansion and genotyping. The method includes the following steps: obtaining third-generation sequencing data and aligning it to a reference genome; extracting STR region sequences based on the coordinate information of the STR in the reference genome and the alignment results; aligning the extracted STR region sequences with a template sequence containing multiple preset STR repeat units and calculating the repeat count; and / or using a method based on adjacent coordinate pairs with the same kmer to calculate the repeat units and repeat counts in the STR region sequences de novo; calculating the repeat count peak using a Gaussian mixture model based on the repeat counts of the STR sequences in multiple reads to obtain the genotyping of the STR in the sample. The method provided by this invention is faster than repeatHMM and also allows for the calculation of the changed repeat unit motif and repeat count when the repeat unit motif changes.

Description

Technical Field

[0001] This invention belongs to the field of gene sequencing technology, specifically relating to detection methods, electronic devices, and storage media for short tandem repeat sequence expansion and genotyping. Background Technology

[0002] Short tandem repeats refer to repetitive sequences in a DNA sequence consisting of multiple nucleotides (with a repeat number greater than or equal to 2 and less than or equal to 6) linked together end to end. Changes in the number of repeat units can have a significant impact on genome structure, which may in turn affect gene expression, modification, and corresponding physiological functions. An increase in the number of short tandem repeat units is called short tandem repeat expansion.

[0003] Second-generation sequencing (NGS) technologies struggle to sequence high-GC repetitive regions and analyze STR repetitive regions with a total length exceeding the read length. Third-generation sequencing (NGS), on the other hand, exhibits no GC bias and can sequence reads longer than most STR repetitive regions, giving it a significant advantage in identifying the number of repetitions in STR expansion. Existing STR analysis tools based on NGS data include Straglr, repeatHMM, tandem-genotypes, and PACMONSTR.

[0004] However, existing STR analysis tools generally have the following drawbacks:

[0005] (1) The process is slow. The three steps of extracting all reads containing STR regions, calling TRF to calculate tandem repeats from scratch, and using BWA-MEM to re-align reads to the genome are slow.

[0006] (2) When the serial repeating region in read expands significantly, the serial repeating region in read cannot be accurately extracted and the number of repetitions cannot be calculated.

[0007] (3) The case where the motif of the repeating unit changes is not considered. Summary of the Invention

[0008] This invention addresses the technical problems existing in the prior art by providing a method for detecting short tandem repeat expansion and genotyping. It is faster than repeatHMM and can calculate the changed motif and the number of repeats when the repeat unit motif changes.

[0009] The present invention adopts the following technical solution to achieve the above technical objectives:

[0010] A method for detecting short tandem repeat expansion and genotyping includes the following steps:

[0011] Step S1: Obtain third-generation sequencing data and use sequence alignment software to align the third-generation sequencing data to the reference genome;

[0012] Step S2: Based on the coordinate information of the STR in the reference genome and the alignment results, extract the corresponding STR region sequence from the read;

[0013] Step S3: Based on the extracted STR region sequence in the read, compare the STR region sequence with a template sequence containing multiple preset STR repeating units, and calculate the number of repetitions.

[0014] and / or

[0015] Using a method based on adjacent coordinate pairs with the same kmer, the repeating units and the number of repeats in the STR region sequence are calculated from scratch, specifically including the following steps:

[0016] Based on all kmers of length K and their coordinates in the extracted STR region sequence of the read, calculate all adjacent coordinate pairs for the same kmer, and determine the minimum sequence between adjacent coordinate pairs as the new repeating unit sequence according to the same rules.

[0017] Based on the difference distribution of all coordinate pairs, filter out coordinate pairs that do not meet the first threshold requirement. Arrange all kmer-filtered coordinate pairs in ascending order. Based on whether the difference between coordinates of the coordinate pairs meets the second threshold and whether the minimum sequence (repeating unit sequence) of the sequence between coordinate pairs is the same, find all main chains that meet the above conditions, where K > 3 and is an integer.

[0018] Based on whether the difference between the start and end coordinates of the main chain meets the third threshold and whether the ratio of the occurrence frequency to the length of the minimum sequence meets the fourth threshold, filter out the main chains that are not within the threshold range, and sort all the filtered main chains in ascending order of coordinates.

[0019] The total number of repetitions is calculated based on the consistency between the repeating unit of the main chain and the minimum sequence of the preset repeating unit. If they are consistent, the repetition count is calculated using the aforementioned comparison method. If they are inconsistent, the total number of repetitions is calculated. The total number of repetitions = the sum of the repetition counts of the main chain + the sum of the sequence repetition counts between each main chain.

[0020] Step S4: Based on the number of repetitions of the STR sequence in multiple reads, the Gaussian mixture model is used to calculate the peak of the number of repetitions, and the genotype of the STR in the sample is obtained.

[0021] This invention can use a method based on adjacent coordinate pairs with the same kmer to calculate the repeating units and the number of repetitions in the STR region sequence from scratch. It is faster than repeatHMM and can also calculate the changed repeating unit motif and the number of repetitions when the repeating unit motif changes.

[0022] Optionally, the method for calculating the number of repetitions by comparing an STR region sequence with a template sequence containing multiple preset STR repeating units also includes adjusting the boundary of the comparison result, removing sequences with poor comparison at both ends of the comparison result, or extending sequences with good comparison to both ends.

[0023] The method for adjusting the alignment result boundaries is as follows: From the alignment results of the STR region sequence and the template sequence, find the start and end coordinate positions (s1, e1) after removing the sequences on both sides. If the number of matching bases in the n bp sequence window at the start or end coordinate is less than the threshold, then shrink inward until the number of matching bases in the n bp sequence window is greater than or equal to the threshold. If the number of matching bases in the n bp sequence window at the start or end coordinate is greater than or equal to the threshold, then extend outward until n If the number of matching bases in the bp window is less than the threshold, the new coordinates (s2, e2) are recorded. If (e2-s2) / (e1-s1) is less than or equal to the set value, the repetition count is calculated to be 0 and returned. Otherwise, the sequence similarity and repetition count of the s2 and e2 intervals are calculated. After removing the gap, the new coordinates (s3, e3) of the s2 and e2 coordinates in the STR region sequence in the read are returned. The repetition count m = c / b, where c is the length of the template sequence s2 containing multiple preset STR repeating units, the length of the sequence without gap in the e2 interval, and b is the length of the repeating unit.

[0024] The number of copies of the repeating units in the template sequence is: or Where a is the length of the STR sequence and b is the length of the repeating unit. If the calculated template sequence length is less than or equal to the STR sequence length, then the number of copies of the repeating unit in the template sequence is: Otherwise, the number of copies of its repeating unit is

[0025] Optionally, the method for determining the minimum sequence is as follows: starting from any base position in the sequence, generate a spliced ​​sequence with the subsequent sequence starting from that base in front of the sequence before that base, and sort them according to the same sorting principle. At the same time, select one of the spliced ​​sequences as the minimum sequence according to the same rule, which is the repeating unit sequence.

[0026] Optionally, the method for determining the main chain includes the following steps: starting from the first coordinate pair among all sorted kmer coordinate pairs, if the difference between the coordinates of the coordinate pairs is within the second threshold range, and the minimum sequence of the coordinate pairs is also the same, then continue searching until the last coordinate pair that meets the condition is found; take the last coordinate pair that meets the condition as the new previous coordinate pair, and continue searching in the aforementioned way until no coordinate pair that meets the condition is found; if no coordinate pair that meets the condition is found, take the range of all coordinate pairs found previously as a main chain, and obtain the main chain (e, f), where e is the starting coordinate of the first coordinate pair, and f is the ending coordinate of the last coordinate pair kmer; start from the next coordinate pair of the main chain and repeat the aforementioned method; if no coordinate pair that meets the condition is found, take the range of all coordinate pairs found previously as a main chain, and repeat the above operation until all main chains are found.

[0027] Optionally, the method for calculating the number of repetitions in the main chain includes the following steps: comparing the main chain sequence with the template sequence that uses the smallest sequence as the repetition unit, calculating the similarity and the number of repetitions; if the number of repetitions is not calculated, then filtering out the main chain.

[0028] The number of repetitions of sequences between main chains is calculated as follows: The repetition units are merged with the expected main chain information and the information of inconsistent main chains and placed in an array, and sorted in ascending order. The ratio of the interval length between main chains to the expected length of the repetition unit is the number of repetitions of sequences between main chains.

[0029] Optionally, the calculation of the total number of repetitions also includes the following steps: adding the length of the interval between the main chains and the length of the two uncovered regions (the remaining read sequences) at the ends of the first and last main chains to calculate the unknown length; if the unknown sequence length divided by the total sequence length is greater than or equal to the fifth threshold, then return 0; otherwise, return the calculated total number of repetitions and record the repetitive units in all main chains.

[0030] Optionally, in step S2, the method for extracting the STR region sequence in the read includes the following steps: extending the first sequence length from the coordinates of the STR region in the reference sequence to both sides to extract the corresponding STR region sequence in the read;

[0031] If the STR region sequence is not extracted, the second sequence length on both sides of the STR region in the reference sequence is extracted and compared to the read. The first similarity threshold is used to determine whether the verification is successful. If the verification is successful, the coordinates of the STR region in the read are recorded. Then, the third sequence length at both ends of the read is extracted and compared with the corresponding region sequence in the reference sequence. The second similarity threshold is used to determine whether the verification is successful. If the verification is successful, the corresponding STR region sequence in the read is extracted. Otherwise, the comparison is skipped.

[0032] The method for detecting short tandem repeat expansion and genotyping provided by this invention, when the expansion of the tandem repeat region in the read is large and the sequence cannot be extracted according to the alignment result in step S1, re-aligns the sequences on both sides of the reference sequence STR with the read sequence, re-extracts the corresponding STR sequence in the read, and then calculates the result of the repeat unit expansion, so that even when the expansion of the tandem repeat region in the read is large, the corresponding STR region in the read can be extracted well, and the result of the repeat unit expansion can be calculated.

[0033] Optionally, in step S2, if the number of STRs is less than the sixth threshold (e.g., 1000), the reads that are aligned with each STR region of the reference sequence are directly extracted from the alignment results of step S1, and then the STR region sequences in the reads are extracted.

[0034] If the number of STRs is greater than or equal to the sixth threshold (e.g., 1000), first sort the alignment results of step S1 according to the reference genome coordinates, sequentially traverse the alignment results, find all STRs within the reference genome range corresponding to each read, and extract the STR region sequence from each STR read.

[0035] Optionally, the STR region sequence is aligned with the template sequence using the Needleman-Wunsch alignment algorithm with a limited band width.

[0036] Using the Needleman-Wunsch algorithm with a band width limit to align the STR region sequence in the read with the template sequence of multiple copies of the repeating unit, and then adjusting the boundary of the alignment result, the STR repetition count information can be calculated more quickly and effectively.

[0037] The present invention also provides an electronic device, including a memory and a processor, wherein the processor is used to execute a computer management program stored in the memory to implement the steps of the above-described method for detecting short tandem repeat expansion and genotyping.

[0038] The present invention also provides a computer-readable storage medium having a computer management program stored thereon, which, when executed by a processor, implements the steps of the methods for detecting short tandem repeat expansion and genotyping as described above. Detailed Implementation

[0039] The specific embodiments of the present invention will be described in further detail below with reference to the examples. These examples are for illustrative purposes only and are not intended to limit the scope of the invention.

[0040] This invention addresses the shortcomings of existing technologies for detecting short tandem repeats. These technologies involve extracting all reads containing STR regions, calculating tandem repeats ab initio using TRF, and employing BWA-MEM for read re-alignment to the genome, which is slow and does not account for changes in repeat unit motifs. This invention provides a novel method for detecting short tandem repeat expansion and genotyping, specifically including the following steps:

[0041] Step S1: Obtain third-generation sequencing data and use sequence alignment software to align the third-generation sequencing data to the reference genome;

[0042] Step S2: Based on the coordinate information of the STR in the reference genome and the alignment results, extract the corresponding STR region sequence from the read;

[0043] Step S3: Based on the extracted STR region sequence in the read, compare the STR region sequence with a template sequence containing multiple preset STR repeating units and calculate the number of repetitions; and / or use a method based on adjacent coordinate pairs with the same kmer to calculate the repeating units and the number of repetitions in the STR region sequence from scratch, specifically including the following steps:

[0044] Based on all kmers of length K and their coordinates in the extracted STR region sequence of the read, calculate all adjacent coordinate pairs for the same kmer, and determine the minimum sequence between adjacent coordinate pairs as the repeating unit sequence according to the same rules.

[0045] Based on the difference distribution of all coordinate pairs, filter out coordinate pairs that do not meet the first threshold requirement. Arrange all kmer-filtered coordinate pairs in ascending order. Based on whether the difference between coordinates of the coordinate pairs meets the second threshold and whether the minimum sequence (repeating unit sequence) of the sequence between coordinate pairs is the same, find all main chains that meet the above conditions, where K > 3 and is an integer.

[0046] Based on whether the difference between the start and end coordinates of the main chain meets the third threshold and whether the ratio of the occurrence frequency to the length of the minimum sequence meets the fourth threshold, filter out the main chains that are not within the threshold range, and sort all the filtered main chains in ascending order of coordinates.

[0047] The total number of repetitions is calculated based on the consistency between the repeating unit of the main chain and the minimum sequence of the preset repeating unit. If they are consistent, the repetition count is calculated using the aforementioned comparison method. If they are inconsistent, the total number of repetitions is calculated. The total number of repetitions = the sum of the repetition counts of the main chain + the sum of the sequence repetition counts between each main chain.

[0048] Step S4: Based on the number of repetitions of the STR sequence in multiple reads, the Gaussian mixture model is used to calculate the peak of the number of repetitions, and the genotype of the STR in the sample is obtained.

[0049] Example 1

[0050] A method for detecting short tandem repeat expansion and genotyping includes the following steps:

[0051] Step S1: Obtain third-generation sequencing data and use sequence alignment software to align the third-generation sequencing data to the reference genome; the input files used include: a BAM format file obtained after aligning the third-generation sequencing reads to the human reference genome using minimap2, STR region coordinate information and repetitive unit sequence information on the provided reference genome, and the reference genome FASTA sequence file.

[0052] Step S2: Based on the coordinate information of the STR in the reference genome and the alignment results, extract the corresponding STR region sequence from the read;

[0053] Specifically, the STR region sequence in the corresponding read is extracted by extending the first sequence length from the coordinates of the STR region of the reference sequence to both sides. The following operation mode can be selected according to the number of STRs: (1) When the number of STRs is less than the sixth threshold, the reads that are aligned with each STR region of the reference sequence are directly extracted from the BAM alignment file, and then the STR region sequence in each read is extracted; (2) If the number of STRs is greater than or equal to the sixth threshold, the BAM format alignment file is sorted according to the coordinates of the reference genome, the read alignment results in the BAM file are sequentially traversed, all STRs in the reference genome range corresponding to each read are found, and the STR region sequence in each read is extracted.

[0054] If the STR region sequence is not extracted, the second sequence length on both sides of the STR region in the reference sequence is extracted and compared to the read. The first similarity threshold is used to determine whether the verification is successful. If the verification is successful, the coordinates of the STR region in the read are recorded. Then, the third sequence length at both ends of the read is extracted and compared with the corresponding region sequence in the reference sequence. The second similarity threshold is used to determine whether the verification is successful. If the verification is successful, the corresponding STR region sequence in the read is extracted. Otherwise, the comparison is skipped.

[0055] When the expansion of the tandem repeat region in the read is large, resulting in the failure to extract the sequence according to the alignment result in step S1, the sequences on both sides of the reference sequence STR are re-aligned with the read sequence, the corresponding STR sequence in the read is extracted again, and then the result of the repeat unit expansion is calculated, so that even when the expansion of the tandem repeat region in the read is large, the result of the repeat unit expansion can be calculated well.

[0056] Step S3: Calculate the number of repetitions based on the extracted STR region sequence from the read;

[0057] S3-1, when motif changes are not considered, the STR region sequence is compared with a template sequence containing multiple preset STR repeating units. The alignment result boundaries are adjusted, and sequences with poor alignment at both ends or sequences with good alignment extending towards both ends are removed. Similarity and repetition count are calculated. If the similarity is lower than the third similarity threshold, it is considered that the repetition count of the STR sequence has not been calculated, and the STR sequence in that read is skipped. The copy count of the repeating units in the template sequence is... or Where a is the STR sequence length and b is the length of the repeating unit, if the calculated template sequence length is less than or equal to the STR sequence length, then the number of copies of the repeating unit in the template sequence is: Otherwise, the number of copies of its repeating unit is

[0058] The STR region sequence is aligned with a template sequence containing multiple predefined STR repeat units using a bandwidth-limited Needleman-Wunsch alignment algorithm. The STR region sequence is the query sequence, and the template sequence containing multiple copies of the STR repeat units is the target sequence. The bandwidth is set to 101. If the bandwidth is less than the length difference between the two sequences plus 1, the bandwidth is set to the length difference between the two sequences plus 1. If the bandwidth is greater than or equal to 10000, an error is reported and the program exits. Different penalty values ​​are applied to the terminal gap, the gap in the query sequence, and the gap in the target sequence. The output is the query sequence with the gap added, the target sequence with the gap added, and a sequence representing the alignment information for each base pair (gap is represented by -, match by |, and mismatch by .).

[0059] The method for adjusting the alignment result boundaries is as follows: From the alignment results of the STR region sequence and the template sequence, find the start and end coordinate positions (s1, e1) after removing the sequences on both sides. If the number of matching bases in the n bp sequence window at the start or end coordinate is less than the threshold, then shrink inward until the number of matching bases in the n bp sequence window is greater than or equal to the threshold. If the number of matching bases in the n bp sequence window at the start or end coordinate is greater than or equal to the threshold, then extend outward until the number of matching bases in the n bp window is less than the threshold, and record the new coordinates (s2, e2). If (e2-s2) / (e1-s1) is less than or equal to the set value, then calculate the repetition count as 0 and return. Otherwise, calculate the sequence similarity and repetition count of the s2 and e2 intervals. After removing the gap, the new coordinates (s3, e3) of the s2 and e2 coordinates in the STR region sequence in the read are returned. Where the repetition count m = c / b, where c is the number of preset STRs. The template sequence s2 of the repeating unit, the length of the sequence without gaps in the e2 interval, and b is the length of the repeating unit.

[0060] Using the Needleman-Wunsch algorithm with a band width limit to align the STR region sequence in the read with the template sequence of multiple copies of the repeating unit, and then adjusting the boundary of the alignment result, the STR repetition count information can be calculated more quickly and effectively.

[0061] S3-2, When considering the case where the motif changes, for the STR sequence in the extracted read, the method based on adjacent coordinate pairs with the same kmer is used to calculate the repeating units in the STR region sequence from scratch, and to see if they are all consistent with the expected repeating unit sequence or the minimum sequence of the expected repeating unit sequence. If they are inconsistent, it is considered that the motif has changed, and the number of repetitions is calculated based on the result of the calculation from scratch. If the number of repetitions is not calculated, the number of repetitions is calculated according to the case where the motif changes are not considered. If they are consistent, it is considered that the motif has not changed, and the number of repetitions is calculated according to the calculation method described above.

[0062] The method for calculating repeating units in the STR sequence from scratch is as follows:

[0063] a. Based on all kmers of length K and their coordinates in the extracted STR region sequence of the read, calculate all adjacent coordinate pairs for the same kmer (for example, take the two adjacent starting coordinates of the kmer as a coordinate pair and calculate all coordinate pairs), and determine the minimum sequence between adjacent coordinate pairs as the new repeating unit sequence according to the same rules.

[0064] The method for determining the minimum sequence is as follows: Starting from any base position in the sequence, generate a concatenated sequence with the subsequent sequence preceding that base and the sequence preceding that base, and sort these concatenated sequences according to the same sorting principle. Then, select one of these concatenated sequences as the minimum sequence according to the same rules. For example, sort all the concatenated sequences from the first character to the last character in the string according to the alphabetical order of their bases, and select the minimum value among all the concatenated sequences as the minimum sequence. It is understandable that other methods can also be used to determine the minimum sequence, as long as the determination method remains consistent.

[0065] b. Filter out coordinate pairs that do not meet the first threshold requirement based on the difference distribution of all coordinate pairs;

[0066] c. Arrange all kmer-filtered coordinate pairs in ascending order. Based on whether the difference between the coordinates of the coordinate pairs meets the second threshold and whether the minimum sequence (i.e. the repeating unit sequence) of the sequences between the coordinate pairs is the same, find all main chains that satisfy the above conditions, where K > 3 and is an integer.

[0067] The method for determining the main chain includes the following steps: Start searching from the first coordinate pair among all sorted kmer coordinate pairs. If the difference between the coordinates of the coordinate pairs is within a threshold range and the minimum sequence of the coordinate pairs is also the same, continue searching until the last coordinate pair that meets the condition is found. Take the last coordinate pair that meets the condition as the new previous coordinate pair and continue searching in the aforementioned way until no coordinate pair that meets the condition is found. If no coordinate pair that meets the condition is found, take the range of all coordinate pairs that have been found as a main chain to obtain the main chain (e, f), where e is the starting coordinate of the first coordinate pair and f is the ending coordinate of the last coordinate pair kmer. Repeat the aforementioned method starting from the next coordinate pair of the main chain. If no coordinate pair that meets the condition is found, take the range of all coordinate pairs that have been found as a main chain and repeat the above operation until all main chains are found.

[0068] d. Based on whether the difference between the starting and ending coordinates of the main chain meets the third threshold and whether the ratio of the occurrence frequency of the minimum sequence to its length meets the fourth threshold, filter out the main chains that are not within the threshold range, and sort all the filtered main chains in ascending order of coordinates.

[0069] e. Align the main chain sequence with the template sequence that uses the minimum sequence as the repeating unit, calculate the similarity and the number of repetitions. If the number of repetitions is not calculated (e.g., the similarity is lower than the third similarity threshold), then filter out the main chain.

[0070] f. If the repeating unit of the main chain is inconsistent with the minimum sequence of the preset repeating unit, then calculate the total number of repetitions. The total number of repetitions = the sum of the number of repetitions of the main chain + the sum of the number of repetitions of the sequences between each main chain; otherwise, use the aforementioned comparison method (method S3-1) to calculate the number of repetitions.

[0071] The method for calculating the total number of repetitions includes the following steps:

[0072] The repetition count of each main chain is calculated as follows: the main chain sequence is compared with a template sequence containing multiple main chain repeating units, and the repetition count is calculated.

[0073] The number of repetitions of sequences between main chains is calculated as follows: The repetition units are merged with the expected main chain information and the information of inconsistent main chains and placed in an array, and sorted in ascending order. The ratio of the interval length between main chains to the expected length of the repetition unit is the number of repetitions of sequences between main chains.

[0074] The calculation of the total number of repetitions also includes the following steps: the interval length between the main chains and the length of the two uncovered areas (the remaining read sequence) at the ends of the first and last main chains are added together to calculate the unknown length. If the unknown sequence length divided by the total sequence length is greater than or equal to the fifth threshold, then 0 is returned; otherwise, the calculated total number of repetitions is returned, and the repetitive units in all main chains are recorded.

[0075] This invention can use a method based on adjacent coordinate pairs with the same kmer to calculate the repeating units and the number of repetitions in the STR region sequence from scratch. This is faster than calling TRF for calculation and also allows the calculation of the changed repeating unit motif and the number of repetitions when the repeating unit motif changes.

[0076] Step S4: Based on the number of repetitions of the STR sequence in multiple reads, a Gaussian mixture model is used to calculate the peak repetition count to obtain the genotype of the STR in the sample. Specifically: Based on the number of repetitions of the STR sequence in multiple reads, a Gaussian mixture model is used with the n_component parameter set to 1 or 2 to fit the distribution of the repetition counts, and the AIC (Akaike information criterion) is calculated. The model with the smaller AIC is selected. If the selected model has an n_component of 1, the average peak value of the model is directly returned as the number of repetitions in the sample. If the selected model has an n_component of 2, the number of repetitions near the average peak value of the model and with a read count greater than or equal to the minimum read count threshold (2) is selected as the number of repetitions of the STR in the sample.

[0077] This embodiment also provides an electronic device, including a memory and a processor, wherein the processor is used to execute a computer management program stored in the memory to implement the steps of the above-described method for detecting short tandem repeat expansion and genotyping.

[0078] This embodiment also provides a computer-readable storage medium storing a computer management program thereon, which, when executed by a processor, implements the steps of the above-described methods for detecting short tandem repeat expansion and genotyping.

[0079] Example 2

[0080] The main objective of the GrandSTR software development in this invention is to detect the number of repeats of STRs in known STR regions of the human genome based on third-generation sequencing technology, and to calculate the genotyping information of STRs in a sample. This embodiment provides a method for detecting short tandem repeat expansion and genotyping, comprising the following steps:

[0081] Step S1: Obtain third-generation sequencing data and use sequence alignment software to align the third-generation sequencing data to the reference genome.

[0082] The input files used in this embodiment include: a BAM format file obtained by aligning third-generation sequencing reads to the human reference genome using minimap2, STR region coordinate information and repeat unit sequence information on the provided reference genome, and a reference genome FASTA sequence file.

[0083] Step S2: Based on the coordinate information of the STR in the reference genome and the alignment results, extract the corresponding STR region sequence from the read.

[0084] In this embodiment, firstly, based on the coordinate information of the STR in the reference genome and the BAM format alignment file obtained after alignment, the STR region sequence in the corresponding read is extracted by extending 30bp to both sides of the STR region coordinate of the reference sequence. Specifically, the following operating mode can be selected according to the number of STRs: (1) When the number of STRs is less than the sixth threshold (e.g., 1000), the reads aligned with each STR region of the reference sequence are directly extracted from the BAM alignment file, and then the STR region sequence in each read is extracted; (2) If the number of STRs is greater than or equal to the sixth threshold (e.g., 1000), the BAM format alignment file is sorted according to the coordinates of the reference genome, the read alignment results in the BAM file are sequentially traversed, all STRs within the reference genome range corresponding to each read are found, and the STR region sequence in each read is extracted.

[0085] If the STR region sequence is not extracted from the BAM format alignment file, the reference sequence and the read sequence are re-aligned. The 1kb sequences flanking the STR region in the reference sequence are extracted and aligned to the read using edlib. The similarity is checked to be greater than 0.75, and the coordinates of the STR region in the read are recorded. If the check is successful, the 1kb sequences at both ends of the read are extracted and aligned with the corresponding 5kb upstream and downstream sequences in the reference sequence using edlib. The similarity is checked to be greater than 0.75. If both checks are successful, the corresponding STR region sequence in the read is extracted. If the checks fail, the alignment is skipped.

[0086] Step S3: Calculate the number of repetitions of repeating units in the extracted STR sequence based on the STR sequence in the read.

[0087] S3-1, When the motif changes are not considered, the number of repetitions of repeating units in the STR sequence is calculated as follows:

[0088] The extracted STR sequence is aligned with a template sequence containing multiple copies of STR repeating units. The alignment boundary is adjusted, removing poorly aligned sequences at either end or well-aligned sequences extending towards either end. Similarity and the number of repetitions are calculated. If the similarity is below the third similarity threshold of 0.75, the STR sequence is considered not to have had its repetition count calculated, and the STR sequence in that read is skipped. The number of copies of repeating units in the template sequence is... or Where a is the STR sequence length and b is the length of the repeating unit, if the calculated template sequence length is less than or equal to the STR sequence length, then the number of copies of the repeating unit in the template sequence is: Otherwise, the number of copies of its repeating unit is

[0089] The alignment of the STR region sequence with the template sequence of multiple copies of the STR repeat unit uses the Needleman-Wunsch alignment algorithm. The STR region sequence is the query sequence, and the template sequence of multiple copies of the STR repeat unit is the target sequence. The band width is set to 101. If the band width is less than the length difference between the two sequences + 1, the band width is set to the length difference between the two sequences + 1. If the band width is greater than or equal to 10000, an error is reported and the program exits. Different penalty scores are used for the terminal gap, the gap in the query sequence, and the gap in the target sequence. The output results are the query sequence with the gap added, the target sequence with the gap added, and the sequence representing the alignment information for each base pair (gap is represented by -, match by |, and mismatch by .).

[0090] The method for adjusting the alignment result boundary is as follows: Input the query and target sequences with gaps calculated by the Needleman-Wunsch alignment algorithm. Find the start and end coordinate positions (s1, e1) of the query and target sequences after removing the sequences on both sides. If the number of matching bases in the 5bp sequence window at the start or end coordinate is less than the threshold (3.95), then shrink inward until the number of matching bases in the 5bp sequence window is greater than or equal to the threshold 3.95. If the number of matching bases in the 5bp sequence window at the start or end coordinate is greater than or equal to the threshold 3.95, then extend outward until the number of matching bases in the 5bp window is less than the threshold, and record the new coordinates (s2, e2). If (e2-s2) / (e1-s1) is less than or equal to the set value 0.3, then the number of repetitions is calculated to be 0, and the result is returned. Otherwise, the sequence similarity and number of repetitions in the S2 and e2 intervals are calculated. The new coordinates (s3, e3) of the S2 and e2 coordinates in the STR region sequence in the read after removing the gap are returned.

[0091] Where, the number of repetitions n = c / b, c is the template sequence s2 of the STR repeating unit copied multiple times, e2 is the length of the sequence without gaps, and b is the length of the repeating unit.

[0092] S3-2, When considering the case where the motif changes, for the STR sequence in the extracted read, the method based on adjacent coordinate pairs with the same kmer is used to calculate the repeating unit sequence from scratch to see if they are consistent with the expected minimum repeating unit sequence. If they are inconsistent, it is considered that the motif has changed, and the number of repetitions is calculated based on the result of the calculation from scratch. If the number of repetitions is not calculated, the number of repetitions is calculated according to the case where the motif changes are not considered. If they are consistent, it is considered that the motif has not changed, and the number of repetitions is calculated according to the calculation method described above.

[0093] The method for calculating repeating units in the STR sequence from scratch is as follows:

[0094] a. Calculate adjacent coordinate pairs and repeating unit sequences. Calculate all kmer of length 5 in the sequence and their coordinates. For the same kmer, calculate all adjacent coordinate pairs. Calculate the minimum sequence based on the sequence between the two coordinates of adjacent coordinate pairs and use it as the new repeating unit sequence. That is, for each base position of the sequence between the two coordinates of a coordinate pair, generate a concatenated sequence with the subsequent sequence starting from that base first and the sequence preceding that base last. Calculate the string values ​​of all concatenated sequences in alphabetical order, and take the concatenated sequence with the minimum value as the minimum sequence.

[0095] For example 1: If the repeating unit is ATGCT, then the splicing sequence of the repeating unit is: ATGCT, TGCTA, GCTAT, CTATG, TATGC, and the minimum sequence is: ATGCT.

[0096] For example, 2: ACTGA ACTGA AACTGA, the sequence between coordinate pairs (0, 5) is ACTGA, the spliced ​​sequences are: ACTGA, CTGAA, TGAAC, GAACT, AACTG, and the smallest sequence is: AACTG.

[0097] b. Filter all adjacent coordinate pairs of kmer. Merge all coordinate pairs of kmer into a list, calculate the distribution of the coordinate differences between the two coordinates of all adjacent coordinate pairs of kmer, take 1.5 times the coordinate difference at the peak of the distribution as the threshold, filter out coordinate pairs with two coordinate differences greater than the threshold, and sort the filtered coordinate pairs in ascending order;

[0098] c. Calculate the main chain based on adjacent coordinate pairs. Initially, the first coordinate pair is taken as the previous coordinate pair. Starting from the coordinates after the previous coordinate pair, if the difference between the coordinates of the current coordinate pair and the previous coordinate pair is within the threshold range (4 times the threshold in step b), and the minimum sequence of the current coordinate pair is consistent with the minimum sequence of the previous coordinate pair (i.e., the first coordinate pair), then continue searching until the last coordinate pair that meets the condition is found. The last coordinate pair that meets the condition is taken as the new previous coordinate pair, and the search continues in the aforementioned way until no coordinate pair that meets the condition is found. If no coordinate pair that meets the condition is found, the range of all coordinate pairs found previously is taken as a main chain. Starting from the next coordinate pair of the main chain, the aforementioned method is repeated to find all main chains.

[0099] The following example, using a specific sequence, illustrates in detail the method for determining the main chain:

[0100] ...CAGCAGCAGCAGCACCAGCAGCAGCTGCAGCAG...

[0101] The included kmer of length 5 and its coordinate pairs are as follows.

[0102] CAGCA:(0,3),(3,6),(6,9),(9,15),(15,18),(18,27), for example, the concatenation sequence of the sequence (CAG) between coordinate pairs (0,3) is CAG, AGC, GCA, sorted in alphabetical order, with the smallest sequence being AGC;

[0103] AGCAG: (1,4), (4,7), (7,10), (16,19), (19,28). For example, the splicing sequences of the sequence (AGC) between (1,4) are AGC, GCA, and CAG. Sorted alphabetically, the smallest sequence is AGC.

[0104] GCAGC:(2,5),(5,8),(8,17),(17,20),(20,26), for example, the splicing sequences of the sequence (GCA) between (2,5) are GCA, CAG, and AGC. Sorted alphabetically, the smallest sequence is AGC.

[0105] The peak difference of all coordinate pairs is 3. The result of merging and sorting all the kmer coordinate pairs after filtering is: (0,3),(1,4),(2,5),(3,6),(4,7),(5,8),(6,9),(7,10),(15,18),(16, 19),(17,20). The threshold is 18, and the minimum sequence is AGC.

[0106] The search begins with the initial coordinate pair (0, 3) as the previous coordinate pair. The last matching coordinate pair is (17, 20), and the new previous coordinate pair is (17, 20). The search continues in the same manner until no matching coordinate pair is found. If no matching coordinate pair is found, the range of all previously found coordinate pairs is taken as a main chain. The main chain is the sequence between coordinate pairs (0, 3) and (17, 20), i.e., the sequence within the range [0, 20]. Considering that the aforementioned coordinates are the starting coordinates of k-mer, and the ending coordinate of a k-mer starting at coordinate 20 is 24, the range of the main chain sequence is [0, 24]. The search continues from the next coordinate pair of the main chain, repeating the above method to find all main chains.

[0107] d. Filter the main chain and calculate the repetition count of repeating units in the main chain. For each main chain, calculate the start and end coordinates, and calculate the number of times the pattern of the smallest sequence appears in the main chain. If the difference between the end coordinate and the start coordinate is less than or equal to the threshold (the threshold in step b), ignore the main chain. If the number of times the pattern of the smallest sequence appears divided by the pattern length is less than the threshold (1.5), also ignore the main chain. Extract the sequence between the start and end coordinates as the query sequence, and the pattern of the aforementioned smallest sequence as the repeating unit. Compare the query sequence with the template sequence according to the method in S3-1 above, adjust the boundary of the comparison result, and calculate the similarity, the number of repetitions, and the adjusted start and end coordinates. If the number of repetitions is not calculated (e.g., the similarity is lower than the third similarity threshold), skip the main chain.

[0108] e. Read the repeating units and main chain information calculated from the beginning. If the repeating unit (i.e., the minimum sequence) is consistent with the minimum sequence of the expected repeating unit (i.e., the preset repeating unit), return 0 (i.e., calculate the number of repetitions according to the method of S3-1). Otherwise, merge the repeating units with the expected and inconsistent main chain information into an array, and sort them according to the start and end coordinates of the main chain information. Add the number of repetitions for each main chain, and add the number of repetitions for the interval between main chains calculated by dividing the interval length by the expected repeating unit length, to obtain the total number of repetitions. Add the length of the interval between main chains and the length of the two uncovered areas at the ends of the first and last main chains (i.e., the remaining read sequences not included in the main chains at both ends) to calculate the unknown length. If the unknown length divided by the total sequence length (the total length of the extracted read sequences) is greater than or equal to the threshold (0.3), return 0 (i.e., calculate the number of repetitions according to the method of S3-1). Otherwise, return the calculated total number of repetitions and record the repeating unit pattern in all main chains.

[0109] The following section provides a detailed explanation of the method for calculating the number of repetitions, using specific examples:

[0110] For example, the extracted read sequence is

[0111]

[0112] The two bolded sequences represent the two main chain sequences, with repetition counts of 4 and 5 respectively. The interval sequence between the two main chains has a length of 6 and a repetition count of 2. Therefore, the total number of repetitions is 4 + 2 + 5 = 11.

[0113] Step S4: Based on the number of repetitions of the STR sequence in multiple reads, a Gaussian mixture model is used to calculate the peak number of repetitions, thus obtaining the genotyping of the STR in the sample. Based on the number of repetitions of the STR sequence in multiple reads, a Gaussian mixture model is used with the n_component parameter set to 1 or 2 to fit the distribution of repetition counts, and the AIC (Akaike information criterion) is calculated. The model with the smaller AIC is selected. If the selected model has an n_component of 1, the average peak value of the model is directly returned as the number of repetitions in the sample. If the selected model has an n_component of 2, the number of repetitions near the average peak value of the model and with a read count greater than or equal to the minimum read count threshold (2) is selected as the number of repetitions of the STR in the sample.

[0114] Example 3

[0115] Information on the coordinates, repeat units, and normal and expanded repeat counts of 47 expanded STRs in the human hg19 reference genome was compiled from the literature. For each STR, one repeat count was randomly selected from both the normal and expanded repeat count ranges. Two sets of STR region sequences were generated based on the repeat unit sequence and the two repeat counts: normal and expanded. Based on the STR's coordinates in the reference genome, 100kb sequences on either side of the coordinates were cropped and added to both sides of the aforementioned two sets of STR region sequences, resulting in two sets of normal and expanded STR reference sequences. Based on the FASTA file of these two sets of STR reference sequences, 15× simulation data was generated using pbsim, equivalent to obtaining a set of 30× simulation data for the reference sequences. The two repeat counts of each STR were recorded in the sium_infos file, based on the pa file containing the coordinates, repeat units, and other information of the aforementioned two sets of STRs in the reference sequences.

[0116] Minimap2 was used to align 30× simulated data to a reference sequence FASTA file, resulting in a BAM file of the alignment results. The reference sequence FASTA file, the BAM file of the alignment results, and a PA file containing coordinates and repeat unit information were used as input files. GrandSTR, as described in Example 1, was used to calculate the number of repeats of each STR in each set of reads, and the genotype of the simulated data samples was calculated. The output STR_infos file contained the STR name, positive and negative strands, the number of STR repeats in the reference sequence read from the PA file, the calculated number of STR repeats in the reference sequence, the STR similarity in the reference sequence, the total number of reads, the genotype, the detailed distribution of repeat counts, the motif indicating a change in repeat unit, the genotype of a changed repeat unit, and the detailed distribution of repeat counts where a change in repeat unit also occurred.

[0117] Repeat the above process 100 times, and calculate the correlation between the actual number of repetitions in the simu_infos file and the genotypes of the simulated data samples. The correlation value is 0.4698. Using the same dataset, the existing repeatHMM, tandem-genotypes, and pacmonSTR methods were used to calculate the genotypes of the simulated data samples, and the correlations with the actual values ​​were 0.0258, 0.2460, and 0.2230, respectively. Specifically, after calculating the number of repetitions for each read using tandem-genotypes, the genotyping statistics code of GrandSTR was used to calculate the genotypes for result comparison.

[0118] The total time (CPU hours) consumed by the four software programs GrandSTR, repeatHMM, tandem-genotypes, and pacmonSTR to run the aforementioned 100 simulation tests were 4.69, 25.90, 87.56, and 246.17, respectively. These results indicate that when the tandem repeat region expansion in the read is large, GrandSTR can better calculate the number of repeats and genotyping results compared to software such as repeatHMM, and it also runs faster.

[0119] Example 4

[0120] Download the HG002 ONT ultralong dataset and align it to the human hg19 reference genome using minimap2 to obtain a BAM file. Extract simple repeat information (3-6 repeat units) from the human repeatmasker results downloaded from the UCSC website to obtain a PA file containing STR genomic coordinates, repeat unit information, etc. Use the hg19 reference genome sequence FASTA file, the alignment result BAM file, and the PA file with coordinates and repeat unit information as input files. Use GrandSTR to calculate the repeat count of each STR in each read and calculate the genotype of the data sample.

[0121] The output STR_infos file contains the STR name, positive and negative strands, the number of STR repeats in the reference sequence read from the pa file, the calculated number of STR repeats in the reference sequence, the STR similarity in the reference sequence, the total number of reads, the genotype, the detailed distribution of repeat counts, the motif of repeated unit alterations, the genotype of repeated unit alterations, and the detailed distribution of repeat counts for repeated unit alterations. A detailed repeat count distribution of NA in the STR_infos file indicates that the repeat count of this STR was not calculated in the reads. This could be due to the absence of a repeat unit sequence for this STR in the corresponding region of the read, or a low similarity between the found sequence and the STR template sequence. If the detailed repeat count distribution for repeated unit alterations in the STR_infos file is not NA, it indicates that the repeat unit of this STR has been altered in the reads, and the repeat count information is derived from de novo calculations.

[0122] Of the 107,287 input STRs, 16,554 STRs have a "Detailed Distribution of Repeats" value of NA, indicating that the repeat counts of these STRs were not calculated. Of the 16,554 STRs, 1,749 STRs have a "Detailed Distribution of Repeats with Changed Repeat Motifs" value that is not NA, indicating that the repeat motifs of these STRs have changed.

[0123] For the same dataset, repeatHMM was used to calculate the number of repetitions of each STR in each read, and the genotypes of the data samples were calculated. Results with a genotype of 0 were removed from the repeatHMM results, and results with a "detailed distribution of repetitions" value of NA were removed from the GrandSTR results. The correlation between the remaining results and the genotype data was 0.8591, indicating that the calculation results of GrandSTR and repeatHMM are quite similar.

[0124] The method for detecting short tandem repeat expansion and genotyping (GrandSTR) provided by this invention uses minimap2 to align reads to a reference genome to obtain a BAM format file. When there are many input STR regions, the BAM file is sorted according to the coordinates of the reference genome, and the BAM file is traversed smoothly to find all STRs corresponding to each read range. The STR sequence in each read is extracted and calculated. The STRs in the sequence are calculated from scratch using a self-developed algorithm. At the same time, the Needleman-Wunsch algorithm with a band width limit is used to compare the STR region sequence in the read with the template sequence of the repeat unit copied multiple times. Then the boundary is adjusted and the number of STR repeats is calculated. The running speed is faster than existing methods.

[0125] When the expansion of the concatenated repeat region in the read is large, resulting in the failure to extract the sequence from the bam file, using edlib to re-align the STR sequences on both sides of the reference sequence with the read sequence and re-extract the corresponding STR sequence from the read can better calculate the result of the repeat expansion.

[0126] For the extracted STR sequence in the read, the repeating units are calculated from scratch to see if they are consistent with the expected repeating units. If they are inconsistent, it is assumed that the motif of the repeating unit has changed. The number of repetitions is calculated based on the result of the calculation from scratch, which fully considers the case where the motif of the repeating unit has changed.

[0127] It should be noted that the descriptions of each embodiment in the above embodiments have different focuses. For parts that are not described in detail in a certain embodiment, please refer to the relevant descriptions in other embodiments.

[0128] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0129] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method of detecting short tandem repeat expansion and genotyping, characterized by, The method comprises the following steps: Step S1, obtaining sequencing data and aligning the sequencing data to a reference genome using sequence alignment software; Step S2, extracting the STR region sequence in the read according to the coordinate information of the STR in the reference genome and the alignment result; Step S3, calculating the number of repetitions according to the extracted STR region sequence in the read; When the motif is not considered to change, the STR region sequence is aligned with a template sequence containing a plurality of preset STR repeat units to calculate the number of repetitions; When the motif is considered to change, The method for calculating the repeat unit and the number of repetitions in the STR region sequence from scratch using the same kmer-based adjacent coordinate pair comprises the following steps: According to all kmers with a length of K and their coordinates in the extracted STR region sequence in the read, all adjacent coordinate pairs of the same kmer are calculated, and the minimum sequence between the adjacent coordinate pairs is determined as a new repeat unit sequence according to the same rule; According to the difference value distribution of all coordinate pairs, coordinate pairs that do not meet the first threshold requirement are filtered out, and the filtered coordinate pairs of all kmers are arranged in ascending order, and the main chain that meets the above conditions is found according to whether the difference value between the coordinates of the coordinate pairs meets the second threshold value and whether the minimum sequence between the coordinate pairs is the same, wherein K>3 and is an integer; According to whether the difference value between the start coordinate and the end coordinate of the main chain meets the third threshold value and whether the ratio of the number of occurrences of the minimum sequence to the length meets the fourth threshold value, the main chain that is not within the threshold value range is filtered out, and all filtered main chains are sorted in ascending order of coordinates; According to the consistency of the repeat unit of the main chain with the minimum sequence of the preset repeat unit, the total number of repetitions is calculated, if consistent, the number of repetitions is calculated using the aforementioned alignment method, if inconsistent, the total number of repetitions is calculated, and the total number of repetitions = the sum of the number of repetitions of the main chains + the sum of the sequence repetition numbers between the main chains; Step S4, calculating the peak value of the number of repetitions using a Gaussian mixture model according to the number of repetitions of the STR sequence in a plurality of reads to obtain the genotyping of the STR in the sample.

2. The method of claim 1, wherein, The method for calculating the number of repetitions by aligning the STR region sequence with the template sequence containing a plurality of preset STR repeat units further comprises the step of adjusting the alignment result boundary to remove the sequences at both ends of the alignment result or extending the better sequences to both ends.

3. The method for detecting short tandem repeat expansion and genotyping according to claim 1, characterized in that, The method for determining the minimum sequence is as follows: starting from any base position in the sequence, a spliced sequence is generated from the subsequent sequence starting from the base and the sequence after the base before the base, and the spliced sequences are sorted according to the same sorting principle, and one of the spliced sequences is selected as the minimum sequence according to the same rule.

4. The method of claim 1, wherein, The determination method of the main chain comprises the following steps: starting from the first coordinate pair in the sorted coordinate pairs of all kmers, if the difference between the coordinates of the coordinate pair is within the second threshold range and the minimum sequences of the coordinate pair are also the same, continue searching until the last coordinate pair meeting the condition is found; taking the last coordinate pair meeting the condition as a new last coordinate pair, continue searching according to the foregoing method until no coordinate pair meeting the condition is found; if no coordinate pair meeting the condition is found, taking the range of all the coordinate pairs found as a main chain to obtain the main chain (e, f), wherein e is the starting coordinate of the first coordinate pair, and f is the ending coordinate of the kmer of the last coordinate pair; starting from the next coordinate pair of the main chain, repeating the above operation until all main chains are found.

5. The method of claim 1, wherein, The calculation method of the number of main chain repetitions comprises the following steps: comparing the main chain sequence with the minimum sequence as a template sequence of a repeating unit to calculate similarity and the number of repetitions, and filtering out the main chain if the number of repetitions is not calculated. The number of sequences between main chains is calculated in the following manner: combining the repeating units and the main chain information of expected consistency and inconsistency in an array and sorting in ascending order, and the ratio of the interval length between main chains to the preset STR repeating unit length is the number of sequences between main chains.

6. The method of claim 1, wherein, In step S2, the method for extracting the STR region sequence in the read comprises the following steps: The STR region sequence in the corresponding read is extracted by extending the STR region coordinates of the reference sequence to both sides by a first sequence length; If the STR region sequence is not extracted, the sequence of the STR region of the reference sequence on both sides by a second sequence length is extracted and compared to the read, whether the verification is successful is judged according to the first similarity threshold, if the verification is successful, the STR region coordinates in the read are recorded, the sequence of the read on both ends by a third sequence length is extracted and compared to the corresponding region sequence in the reference sequence, whether the verification is successful is judged according to the second similarity threshold, if the verification is successful, the corresponding STR region sequence in the read is extracted, otherwise, the comparison is skipped.

7. The method of claim 1, wherein, In step S2, if the number of STRs is less than the sixth threshold, the reads compared with each STR region of the reference sequence are directly extracted from the comparison result of step S1, and then the STR region sequence in the read is extracted; If the number of STRs is greater than or equal to the sixth threshold, the comparison result of step S1 is sorted according to the reference genome coordinates, the comparison result is sequentially traversed, all STRs in the range of the reference genome corresponding to each read are found, and the STR region sequence in the read is extracted for each STR.

8. The method of claim 1, wherein, The comparison between the STR region sequence and the template sequence is performed by using the needleman-wunsch comparison algorithm with a limited band width.

9. An electronic device, comprising: The device comprises a memory and a processor, and the processor is used to execute the computer management program stored in the memory to realize the steps of the detection method and the genotyping method of short tandem repeat expansion according to any one of claims 1-8.

10. A computer-readable storage medium, characterized in that, A computer program product comprising a computer readable medium having stored thereon a computer management program which, when executed by a processor, implements the steps of the method for detecting short tandem repeat expansion and genotyping according to any one of claims 1 to 8.

Images