Method for representing and generating variation path graph of high-throughput sequencing data

Abstract

The present invention discloses a method for representing a variation path graph of high-throughput sequencing data and a method for generating the same. The variation path graph includes a path segment and a cut point. The path segment is composed of a group of path units connected sequentially from left to right. The path segment represents an arbitrary sequence with a length greater than or equal to 1. The beginning and end of the path segment and the path unit are marked by a cut point. The present invention defines a data structure of a variation path graph, including a path segment, a path unit, and a cut point; and a method for representing the original sequencing data read segment with a variation path graph, proposes a method for identifying a cut point, and a method for establishing a corresponding relationship between a cut point "on the variation path graph" and a "position on the original sequencing data read segment", proposes a method for tracing the cut point, splitting the read segment into small segments covering different path units, and then aggregating them into path units, connecting them into path segments, and finally obtaining a splicing method for the variation path graph.

Description

Technical Field

[0001] The present invention relates to the analysis and application field of high-throughput sequencing data, and in particular to a method for representing a variation path graph of high-throughput sequencing data and a method for generating the same. Background Art

[0002] High-throughput sequencing is a common technology in modern biomedical research and application fields. The result of high-throughput sequencing is a large amount of short fragment sequence information, called "reads". The current common read form is a pair of strings with a length of 100-150. One of the strings consists of A, T, C, and G, representing the measured sequence; the other string consists of characters converted from the sequencing quality values ​​corresponding to each position of the sequence. The analysis process of high-throughput sequencing data generally requires first aligning the sequences of these reads to the corresponding positions of the reference genome. Then, based on how many reads cover a specific position on the reference genome (i.e., coverage depth) and the difference between the reference genome and the measured sequence at the corresponding position, a series of biomedical indicators are calculated, such as: gene expression, transcription factor binding position, epigenetic modification position, gene variation, protein variation, possible new antigens (peptide variation), etc.

[0003] At present, this mainstream analysis process has certain defects, mainly: 1) Due to individual differences in the genome of the sequencing object, some of the genome sequences unique to the sequencing object do not exist in the reference genome, so this part of the sequencing data cannot be matched to the reference genome. As a result, the sequence features related to the unique genome sequences of the sequencing object cannot enter the field of view of the subsequent analysis process, becoming a blind spot in biomedical research and application. This defect is very important in tumor-related research and applications. The cause of tumors is likely to be due to mutations in the individual-specific genomic regions of the patient, and these mutations cannot be discovered using the analysis process based on the reference genome. 2) At present, software tools that match sequencing reads to the reference genome generally have certain design limitations: it is difficult to effectively discover chromosome structural variations by comparing the reference genome method (for example, the genome of the sequencing object has undergone large-scale sequence rearrangements relative to the reference genome, etc.). Chromosome structural variations often have a greater impact on gene function. The inefficiency of the reference genome-based method in discovering chromosome structural variations is also one of the bottlenecks in the current research on tumor and other disease mechanisms.

[0004] On the other hand, the genome of the sequencing object can be spliced ​​through high-throughput sequencing data. This method can obtain the unique genome sequence of the sequencing object. However, due to the complex structure of the genome, existing sequence splicing tools usually only give a set of the most likely gene sequences. These sequences are selected based on empirical possibilities and cannot fully represent all the sequence variation information contained in the real sequencing data. Therefore, the personalized genome sequence of the sequencing object obtained using the current splicing tools still has information loss. In addition, the use of the spliced ​​personalized genome as a reference genome for subsequent analysis is still subject to the second limitation of the above-mentioned reference genome-based analysis process.

[0005] At present, the basic algorithm of the commonly used genome splicing tool is the de Bruijn Graph. In the de Bruijn graph, the sequence string of each read is divided into small segments of length K, called K strings. K strings with the same sequence are considered to be the same locus on the genome (therefore, K is also referred to as the resolution parameter in the present invention). Each K string is a node. If two K strings have overlapping (K-1) characters, there is an edge between them. In this way, each read can be mapped to a path in the de Bruijn graph. Thus, the sequence splicing problem is transformed into the problem of finding a set of paths containing all reads in the de Bruijn graph. Although the de Bruijn graph can be used to fully represent high-throughput sequencing data, the number of nodes of a de Bruijn graph with a sufficiently high resolution (K is large enough) is too large (almost equal to the number of sites in the entire genome), and it is very difficult to traverse the graph; at the same time, the sequence recorded in each node in the graph has K-1 characters repeated with the sequence recorded in its adjacent node, resulting in inefficient use of storage space. Due to the existence of these two problems, direct analysis of sequencing data based on the de Bruijn graph without relying on a reference genome is very inefficient, computationally expensive, and difficult to use in practice.

[0006] The standard format for high-throughput sequencing data is FASTQ. A unit of a FASTQ file represents a read. Each unit takes up 4 lines of space. The first line is the name of the read, the second line is the sequence information of the read, the third line is a "+" sign, and the fourth line is the quality information of the read. The second line has the same length as the fourth line. The sequencing quality of the character at a specific position in the second line corresponds to the character at the same position in the fourth line. For samples with double-end sequencing, one sample has two FASTQ files, in which the reads correspond to the sequences of the head and tail of the detected nucleic acid molecules and their sequencing quality, respectively. Summary of the invention

[0007] In order to overcome the shortcomings of the above-mentioned existing high-throughput sequencing data analysis methods, the present invention designs a variation path diagram representation method and a generation method of high-throughput sequencing data, which can represent the position correspondence between the reads of all high-throughput sequencing data from multiple samples (rather than just the position correspondence between the reads that can be matched to the reference genome) using the variation path diagram. The variation path diagram records the sequence of each "path fragment" in the diagram, which fragments of which reads of which samples are covered. This information can be easily and intuitively docked with existing sequencing data analysis tools for subsequent biological feature analysis. Overall, the variation path diagram is another representation of the original sequencing data, and the information is not lost during the conversion process from the original sequencing data to the variation path diagram. Converting the original sequencing data into a variation path diagram can easily dock with subsequent biological feature analysis, including analysis of gene expression, transcription factor binding position, epigenetic modification position, gene variation, protein variation, possible new antigens (peptide variation), etc.

[0008] The present invention is achieved through the following technical solutions:

[0009] The present invention discloses a specific representation method of a variation path graph of sequencing data, wherein the variation path graph comprises a path segment and a tangent point, wherein the path segment is composed of a group of path units connected sequentially from left to right, the path segment represents an arbitrary sequence with a length greater than or equal to 1, and the beginning and end of the path segment and the path unit are marked by a tangent point.

[0010] As a further improvement, the path segment of the present invention carries a set of additional information: which segment of the sequence of the path segment is covered by which read segment in the original sequencing data, and the correspondence between the segment of the sequence of the path segment and the segment of the sequence of the read segment covering it;

[0011] A path unit represents an arbitrary sequence with a length greater than or equal to 1. It is the basic data structure used to assemble path fragments. It is only used in the intermediate process of calculation and is not presented in the final mutation path graph.

[0012] A cut point is a number that represents the connection between path units. If the tail cut point of a path unit has the same number as the head cut point of another path unit, the two path units are connected from left to right.

[0013] The variation path graph has a resolution parameter K, that is, in the variation path graph, the sequence composed of any continuous (K-1) characters (which can span path fragments) is inconsistent; if in the original sequencing data, two read segments have (K-1) consecutive characters in common, then the two (K-1) length fragments in the two read segments cover the same path of length (K-1) characters in the variation path graph, which may span multiple interconnected path fragments.

[0014] As a further improvement, when the head or tail of a path segment (or path unit) can connect at least two path segments (or path units) or cannot connect any path segment (or path unit), the head or tail tangent point of the path segment (or path unit) is called a bifurcation point of a path.

[0015] The present invention also discloses a specific method for generating a variation path graph of sequencing data, which utilizes the identification and tracking of cut points to establish an association between "the bifurcation point of the path segment (or path unit) in the variation path graph", "the position inside the path segment sequence corresponding to the starting point or end point of the original sequencing read segment (i.e., the starting point or end point of the path unit)" and "the position coordinates of the sequence within the original sequencing read segment", and utilizes this association to achieve efficient distributed splicing.

[0016] As a further improvement, the steps are as follows: first, scan the read sequences of all original sequencing data to generate all possible subsequences of length (K-1), and record the characters on the left and right sides of the subsequences. Then, cluster these subsequences of length (K-1). According to the possibility of the characters on the left or right sides of each type of subsequence of length (K-1), identify the position of the cut point in the original read sequence, split the original read sequence into small fragments marked by the cut point at the head and tail positions, and then group these small fragments according to the head and tail cut points and sequences to generate path units. Then, connect the non-branched path units head to tail to form path fragments, and obtain a variation path graph.

[0017] As a further improvement, the present invention specifically comprises the following steps:

[0018] 1) For each original sequencing read, scan the (K-1) length strings in its sequence from left to right to generate (K-1) string objects;

[0019] 2) Cluster all (K-1) string objects generated in step 1) according to their (K-1) length strings, assign a unique cluster ID (positive integer) to each cluster of (K-1) string objects, and then analyze each cluster of (K-1) string objects to obtain the cut point-read segment position object;

[0020] 3) clustering all the cut point-read position objects generated in step 2) according to their source read numbers; then analyzing the cut point-read position objects of each cluster to generate path unit-read position objects;

[0021] 4) clustering all the path unit-read position objects generated in step 3) according to the triplets of the start cut point ID, the end cut point ID, and the path unit sequence; and then mapping each cluster of path unit-read position objects into a path unit object;

[0022] 5) For all the path unit objects generated in the previous step, perform the operation of merging equivalent tangent point IDs;

[0023] 6) Perform left tangent point alignment operation and right tangent point alignment operation;

[0024] 7) Repeat step 6) until during the execution of step 6), each cluster contains only one path unit object;

[0025] 8) Based on the path unit object generated in step 7), a set of internal tangent points of the path segment is generated, where the internal tangent points of the path segment are connected to only one path unit (i.e., no bifurcation) both forward and backward;

[0026] 9) Mapping the path unit object generated in step 7) into an unconnected path segment object;

[0027] 10) Connect the unconnected path segment objects generated in step 9);

[0028] 11) The fully connected set of path fragment objects generated by the above steps constitutes a variation path graph representation of the original sequencing data; the path fragment object contains the cut point ID information on the left and right sides, as well as which part of the path fragment sequence is covered by which fragment of the original sequencing data read sequence.

[0029] The present invention also discloses another specific method for generating a variation path graph of sequencing data. When the sequencing data set is very large, the massive original sequencing data is divided into different parts, and a variation path graph is constructed for each part of the original sequencing data respectively. The obtained variation path graphs are then merged to obtain a variation path graph of all sequencing data or specific sequencing data.

[0030] As a further improvement, the specific steps of merging two variation path graphs described in the present invention are as follows:

[0031] 1) Give each position of the path segment in the two variant path graphs a unique position number;

[0032] 2) Treat the two variant path graphs as a read segment respectively, scan all the strings of length (K-1) in the two read segments, and generate (K-1) string objects. The string may start from any position in any path segment and cross the boundary of the path segment with any path; the (K-1) string object no longer records the full sequence information of its source read segment, but records the end position of the (K-1) string on the read segment;

[0033] 3) Cluster the (K-1) string objects generated in step 2) according to their (K-1) length character strings, assign a unique cluster ID (positive integer) to each cluster of (K-1) string objects, and then analyze each cluster of (K-1) string objects respectively, and directly use the end position information of the (K-1) string on the source read segment recorded by the (K-1) string object, instead of calculating the end position based on its starting position, to obtain the cut-point-read segment position object; at the same time, the cut-point-read segment position object needs to record the characters on both sides of the cut-point, but does not record the full sequence of the source read segment. The characters on both sides of the cut-point can be obtained based on the information recorded by the (K-1) string object;

[0034] 4) clustering all the cut point-read position objects generated in step 3) according to their source read numbers; then analyzing the cut point-read position objects of each cluster respectively, performing operations to generate path unit-read position objects; when generating the path unit-read position objects, the path unit sequence is not recorded, but the first and last characters of the path unit are recorded;

[0035] 5) clustering all the path unit-read position objects generated in step 4) according to the four-tuple of the starting cut point ID, the ending cut point ID, and the first and last characters of the path unit; then mapping each cluster of path unit-read position objects into a path unit object; when generating the path unit object, the path unit sequence is not recorded, but the first and last characters of the path unit are recorded;

[0036] 6) For all the path unit objects generated in the previous step, perform the operation of merging equivalent tangent point IDs;

[0037] 7) Based on the path unit object generated in step 6), a set of internal tangent points of the path segment is generated, where the internal tangent points of the path segment are tangent points that are connected to only one path unit (i.e., not bifurcated) both forward and backward;

[0038] 8) Mapping the path unit object generated in step 6) to an unconnected path segment object, where the length of the path unit sequence of the path unit can be calculated based on the distance between the source read start position and the source read end position of any record in its source read segment position table;

[0039] 9) Connect the unconnected path segment objects generated in step 8);

[0040] 10) The starting and ending positions of the fully connected path segment object set generated by the above steps must correspond to the starting position or the ending position of a path segment in the original variation path graph, that is, the path segment of the synthesized variation path graph must be the split of the path segment in the original variation path graph. Based on this correspondence, the source read segment mapping table of the path segment in the original variation path graph can be converted into the source read segment mapping table of the path segment in the synthesized variation path graph;

[0041] 11) The set of path fragment objects generated by the above steps is the result of merging the two variation path graphs, forming a variation path graph representation of the merged original sequencing data; the path fragment object contains the cut point ID information on its left and right sides, as well as which part of the path fragment sequence is covered by which fragment of the merged original sequencing data read sequence.

[0042] The present invention has the following innovative features:

[0043] 1) The data structure of the variation path graph is defined, including path segments, path units, and cut points; and the method of using the variation path graph to represent the raw sequencing data reads.

[0044] 2) A method for identifying the cut point and a method for establishing the correspondence between the "position of the cut point on the variation path diagram" and the "position on the original sequencing data read segment" are proposed.

[0045] 3) A method for splicing the variant pathway graph is proposed by tracking the cut points, splitting the read segments into small fragments covering different pathway units, aggregating them into pathway units, connecting them into pathway fragments, and finally obtaining the variant pathway graph.

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

[0047] 1) Since the process of converting raw sequencing data to variation pathways does not require matching the reference genome, the information of the part of the raw sequencing data that cannot be matched to the reference genome will not be lost. Subsequent data analysis based on variation pathways can not only observe sequence features that can be located on the reference genome, but also sequence features that cannot be located on the reference genome (specific to individual samples or groups).

[0048] 2) The variation path diagram can represent the correspondence between reads of sequencing data from multiple groups and samples. Based on which reads of which samples and which fragments cover the sequence of each "pathway fragment", it is possible to compare and report the differences in sequence features between different samples or groups without relying on the coordinate system provided by the reference genome, including sequence mutations, transcript (or peptide) expression differences, epigenetic modification position differences, transcription factor binding site differences, etc.

[0049] 3) Compared with sequencing data in FASTQ format or de Bruijn graph format, the variation path diagram has almost no data redundancy and can directly calculate: a) how many reads of which samples cover a specific sequence, b) the sequencing quality characteristics of these regions of these reads, c) whether there are bifurcations before and after a specific sequence, and d) how many reads from which samples support these bifurcations. The above four characteristics are the analysis basis of almost all existing sequencing data biological feature analysis tools. The variation path diagram provides sufficient basis for subsequent analysis in a compact and intuitive way, which can be easily connected to existing analysis tools to complete subsequent biological feature analysis.

[0050] 4) The method for generating the mutation path graph can be designed as an algorithm that conforms to the Map-Reduce mode (see the specific implementation section for details), so it can be implemented in a concise manner using existing big data analysis tools (such as Spark, etc.). This implementation method has high computational efficiency and can process massive sequencing data that traditional implementation methods cannot process, and can use a larger resolution parameter K for calculation. With the iterative optimization of the underlying big data analysis framework, the execution efficiency of the algorithm designed by the present invention will also continue to improve. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 It is a relationship diagram between the variation path diagram described in the present invention and the original high-throughput sequencing data read segments. DETAILED DESCRIPTION

[0052] The present invention discloses a method for representing a variation path graph of high-throughput sequencing data and a method for generating the same. Figure 1It is a relationship diagram between the variant path graph of the present invention and the original high-throughput sequencing data read segment. The variant path graph includes path segments and cut points. The path segments include path units. The beginning and end of the path segments and the path units are marked by cut points. If the end cut point of a path segment is the same as the first cut point of another path segment, they are connected from left to right. The beginning and end of a path segment can be connected to multiple path segments or no path segment is connected, which is called "path bifurcation". The path segment is connected by "path units". The beginning and end of the path unit are also marked by cut points. If the end cut point of a path unit is the same as the first cut point of another path unit, they are connected from left to right; the path units inside the path segment are connected one-to-one in order without bifurcation. The first and last path units of a path segment may not be connected to other path units, or may be connected to multiple path units belonging to other path segments. A read segment in the original sequencing data covers a unique path in the variation path graph. The starting and ending positions of each read segment must correspond to a tangent point. If there is an identical subsequence of length (K-1) in two read segments, but the left (or right) side characters adjacent to the identical subsequence are different, then the left (or right) side positions of the identical subsequence of the two read segments correspond to the same tangent point in the variation path graph (i.e., the tangent point at the bifurcation of the path, such as Figure 1 If there is an identical subsequence of length (K-1) in two reads, then the identical subsequence covers the same and unique path in the variation path graph.

[0053] According to the above definition, the cut point can be identified by the following method: first scan the read sequence of all the original sequencing data, generate all possible subsequences of length (K-1), and record the characters on the left and right sides of the subsequence, and then cluster these subsequences of length (K-1). If there is no character on the left side of a subsequence, or there are more than 2 different characters, or the start of a read is this subsequence, then there is a cut point on the left side of the subsequence, and each read containing this subsequence has an identical cut point on the left side of the subsequence. Similarly, if there is no character on the right side of a subsequence, or there are more than 2 different characters, or the end of a read is this subsequence, then there is a cut point on the right side of the subsequence, and each read containing this subsequence has an identical cut point on the right side of the subsequence.

[0054] Since the positions of the variation path graph corresponding to the head and tail of the read sequence must be the tangent points, the part between any pair of adjacent tangent points in the variation path graph (a path unit) must be completely covered by a certain read; and the sequence of any read segment must be spliced ​​by a group of path units connected head to tail, and this group of path units may belong to different path segments, such as Figure 1Read 2, read 4 in. Based on the above characteristics, the sequence fragments in the original reads with the same head and tail cut points and the same intermediate sequence must cover the same path unit in the variation path graph. Therefore, the original read sequence can be split into small fragments marked by the cut points according to the position of the cut points in the original read sequence obtained in the previous step; these small fragments are then grouped according to the head and tail cut points and sequences to generate path units, each of which contains information about which fragments of which original reads cover the path unit; then the unbranched path units are connected head to tail to form path fragments (at the same time, information about which parts of the path fragments are covered by which fragments of which original reads is generated), and the variation path graph is obtained.

[0055] The technical solution of the present invention is further described below through specific implementation cases:

[0056] The present invention discloses a specific representation method of a variation path graph of sequencing data, wherein the variation path graph includes path segments and cut points:

[0057] 1) A path segment represents an arbitrary sequence with a length greater than or equal to 1. Each path segment carries a set of additional information: which segment of the sequence of the path segment is covered by which read in the original sequencing data, and the correspondence between the segment of the sequence of the path segment and the segment of the sequence of the read that covers it;

[0058] 2) A tangent point is a number that represents the connection method between path segments. Each path segment has a tangent point at the beginning and end of the sequence. When the ending tangent point number of a path segment (A) is the same as the starting tangent point number of another path segment (B), it means that in the variation path graph, path segment A connects path segment B from left to right. A path segment can be connected to multiple other path segments at its beginning and end (i.e., path bifurcation). The beginning and end of a path segment may not be connected to any other path segments (i.e., the starting and end points of the variation path graph, also called path bifurcation).

[0059] 3) Each path segment consists of a group of "path units" connected sequentially from left to right. Path units are similar to path segments. Each path unit represents an arbitrary sequence with a length greater than or equal to 1, with a tangent point on the left and right. If the ending tangent point number of a path unit (X) is the same as the starting tangent point number of another path unit (Y), it means that in the variation path graph, the X path unit connects the Y path unit from left to right. The starting tangent point of a path segment is the starting tangent point of the leftmost path unit that constitutes the path segment; the ending tangent point of a path segment is the ending tangent point of the rightmost path unit that constitutes the path segment. Each path unit can only uniquely belong to one path segment; within a path segment, all path units are connected sequentially only from left to right, and there is no bifurcation. The path unit is the basic data structure used to assemble path segments. It is only used in the intermediate process of calculation and is not presented in the final variation path graph.

[0060] 4) The variation path graph has a resolution parameter K (K is a natural number), that is, in the variation path graph, any sequence consisting of (K-1) consecutive characters (which can span path segments) is inconsistent; if in the original sequencing data, two read segments have (K-1) consecutive characters in common, then the two (K-1) length segments in the two read segments cover the same path of (K-1) characters in length in the variation path graph. This path may span multiple interconnected path segments. Therefore, the splicing method of the original sequencing data represented by a variation path graph with a resolution parameter of K is equivalent to a de Bruijn graph splicing of the original sequencing data based on a "K string".

[0061] The present invention also discloses a specific method for generating a variation path graph of sequencing data, which uses the tracking of "cut points" to establish an association between "the bifurcation point of the path segment (or path unit) in the variation path graph", "the position of the path segment sequence corresponding to the starting point or end point of the original sequencing read segment (i.e., the starting point or end point of the path unit) in the variation path graph" and "the position coordinates of the sequence in the original sequencing read segment", and uses this association to achieve efficient distributed splicing, which includes the following steps:

[0062] 1) For each raw sequencing read, scan the (K-1) length strings in its sequence from left to right to generate (K-1) string objects. Each (K-1) string object contains the following information:

[0063] a. The scanned sequence string is of length (K-1);

[0064] b. The character on the left side of the string (if the string is the first (K-1) string of its source segment, the character on the left side is recorded as '|');

[0065] c. The characters on the right side of the string (if the string is the last (K-1) string of its source segment, the characters on the right side are recorded as '|');

[0066] d. The source read segment number of the (K-1) string;

[0067] e. The starting position coordinates of the (K-1) string on the source read segment;

[0068] f. When the (K-1) string is the first (K-1) string of its source read segment, all sequence information of its source read segment is recorded, otherwise it is not recorded.

[0069] 2) Cluster all (K-1) string objects generated in step 1) according to their (K-1) length strings, and assign a unique cluster ID (positive integer) to each cluster of (K-1) string objects; then perform the following operations on each cluster of (K-1) string objects to obtain the “cut point-read segment position” object:

[0070] a. If the characters on the left side of the cluster (K-1) string object contain two or more non-'|' characters, or contain a '|' character, then there is a cut point on the left side of the sequence of the cluster (K-1) string object (when it contains two or more non-'|' characters or only contains a '|' character, the cut point corresponds to the bifurcation point in the variation path graph; when it contains a '|' character and a non-'|' character at the same time, the cut point corresponds to the position inside the path fragment sequence in the variation path graph corresponding to the start of the original sequencing read segment, and is also the starting position of a path unit).

[0071] b. If there is a cut point on the left side of the sequence of the (K-1) string object in the cluster, then each (K-1) string object in the cluster is mapped into a cut point-read position object. Each cut point-read position object contains the following information:

[0072] Cut point ID (equal to the negative value of cluster ID) (negative value indicates left cut point);

[0073] Source read segment number (equal to the source read segment number of the (K-1) string object);

[0074] Source read segment position (equal to the starting position of the (K-1) string object on the source read segment);

[0075] The full sequence of the source read segment (recorded when the (K-1) string object records the full sequence of the source read segment, otherwise not recorded).

[0076] c. If the right side characters of the cluster (K-1) string object contain two or more non-'|' characters, or contain a '|' character, then there is a tangent point on the right side of the sequence of the cluster (K-1) string object (when it contains two or more non-'|' characters or only contains a '|' character, the tangent point corresponds to the bifurcation point in the variation path graph; when it contains a '|' character and one non-'|' character, the tangent point corresponds to the end position of the original sequencing read segment inside the path fragment sequence in the variation path graph, and is also the end position of a path unit).

[0077] d. If there is a cut point on the right side of the sequence of the (K-1) string object in the cluster, then each (K-1) string object in the cluster is mapped into a cut point-read position object. Each cut point-read position object contains the following information:

[0078] Cut point ID (equal to cluster ID) (positive value indicates right cut point).

[0079] Source read number (equal to the source read number of the (K-1) string object).

[0080] Source read segment position (equal to the end position of the (K-1) string object on the source read segment, which can be calculated by the start position of the (K-1) string object on the source read segment).

[0081] The full sequence of the source read segment (recorded when the (K-1) string object records the full sequence of the source read segment, otherwise not recorded).

[0082] e. If there is no cut point on the left and right sides of the sequence of the cluster (K-1) string object, the cluster (K-1) string object does not generate any "cut point-read segment position" object.

[0083] 3) Clustering all the cut point-read position objects generated in step 2) according to their source read numbers; then performing the following operations on the cut point-read position objects in each cluster to generate a "path unit-read position" object:

[0084] a. Sort the cluster cut point-read position objects by their source read positions;

[0085] b. From steps 1)-2), it can be obtained that the first (K-1) string object generated by a certain read segment must carry the full sequence information of the read segment; the left side of the (K-1) string object must be a cut point; therefore, the cut point-read segment position object generated by the (K-1) string object must carry the full sequence information of the read segment and its source read segment position is the smallest (ranked first), therefore, the first object of the sorted cut point-read segment position object set must contain the full sequence information of the read segment, and the full sequence of the read segment is recorded as S.

[0086] c. Scan the sorted tangent point-read position objects from front to back, taking a pair of adjacent tangent point-read position objects each time, and record the object at the front as A and the object at the back as B, and map them into a path unit-read position object, which contains the following information:

[0087] Starting cut point ID (equal to the cut point ID of A);

[0088] End cut point ID (equal to B's cut point ID);

[0089] Source read number (equal to the common source read number of the cluster objects);

[0090] Source read start position (equal to source read position of A);

[0091] The source read segment end position (equal to the source read segment position of B);

[0092] Path unit sequence (equal to a substring extracted from S with the "source read position" of A and B as the starting and ending positions).

[0093] 4) Cluster all the path unit-read position objects generated in step 3) according to their (start cut point ID, end cut point ID, path unit sequence) triples; then map each cluster of path unit-read position objects into a "path unit" object, which contains the following information:

[0094] a. Starting cut point ID (equal to the common starting cut point ID of the cluster objects);

[0095] b. End cut point ID (equal to the common end cut point ID of the cluster objects);

[0096] c. Path unit sequence (equal to the common path unit sequence of clustered objects);

[0097] d. Source read position table. Each path unit-read position object in each cluster is mapped to a record in this table. Each record contains the following information:

[0098] Source read number (equal to the source read number of the path unit-read position object);

[0099] Source read start position (equal to the source read start position of the path unit-read position object);

[0100] Source read end position (equal to the source read end position of the path unit-read position object);

[0101] 5) For all the path unit objects generated in the previous step, perform the operation of merging equivalent tangent point IDs. The steps are as follows:

[0102] a. For all path unit objects, according to the length of their path unit sequence, split them into an object set with non-zero length (recorded as a non-empty path unit set) and an object set with zero length (recorded as an empty path unit set). If the empty path unit set is an empty set, there is no need to merge them, and all path unit objects are directly returned as merged path unit objects;

[0103] b. For the path units in the empty path unit set, the starting point cut ID is equivalent to the ending point cut ID. If two point cut IDs are equivalent to the third point cut ID at the same time, then the three point cut IDs are equivalent to each other. According to this equivalence rule, a set of "equivalent point cut ID sets" can be calculated, and the point cut ID replacement rule is generated: replace the point cut ID in each "equivalent point cut ID set" with the point cut ID with the smallest value in the set; the point cut ID that is not in the set of equivalent point cut IDs remains unchanged;

[0104] c. For the path unit objects in the non-empty path unit set, apply the above-mentioned point-cut ID replacement rule to replace its starting point-cut ID and ending point-cut ID;

[0105] d. Cluster the replaced non-empty path unit objects according to the triples of (starting cut point ID, ending cut point ID, path unit sequence); then map each cluster of path unit objects to a path unit object. The mapping method is: use any path unit object in the cluster as a template for the mapping result, change its source read position table to the union of the records of the source read position tables of all path units in the cluster, and keep the rest of the information unchanged;

[0106] e. The path unit object obtained in the above steps is the merged path unit object.

[0107] 6) For all the path unit objects generated in the previous step, perform the left tangent point alignment operation. From steps 1)-2), it can be obtained that the path bifurcation point in the final variation path graph must correspond to a tangent point. Therefore, the first character on the right side of a tangent point determines the unique path between two adjacent tangent points in the variation path graph. Based on this principle, the path unit can be left-aligned. The left alignment of the path unit can split the fused path unit. The existence of the fused path unit is due to the truncation effect on the left and right sides of the read segment, which causes some tangent points within the length (K-1) on the left and right sides of the read segment to be unable to be marked by steps 1)-2), so that the path unit obtained in steps 3)-4) is not the smallest path unit, but a fusion of several path units that cross the unmarked tangent points. The left alignment operation steps of the path unit are as follows:

[0108] a. For all path unit objects, cluster them according to their (starting point ID, the first character of the path unit sequence) tuples;

[0109] b. Map each cluster of path unit objects to a new cluster of path unit objects as follows:

[0110] 1. Sort the path unit objects in the cluster according to the length of the path unit sequence. The path unit object with the smallest length is A, and its length is L.

[0111] 2. Map A to a new path unit object. The information of the new object is consistent with A, and only the source read segment position table of the object is changed and updated. The source read segment position table of the new object contains:

[0112] i) All records in the source read position table of A;

[0113] ii) Records mapped from records in the source read position table of other path units in the cluster, the mapping rule is: the source read segment number remains unchanged, the source read segment start position remains unchanged, and the source read segment end position is set to a position L (backward) away from the source read segment start position.

[0114] 3. Map the non-A path unit object in the cluster to a new path unit object, and the starting cut point ID of the new object is set to the ending cut point ID of A. The ending cut point ID of the new object remains unchanged. The path unit sequence of the new object is set to the path unit sequence of the original object after deleting the path unit sequence of A from the left (the sequence may be empty). The records of the source read position table of the new object are mapped from the records of the source read position table of the original object. The mapping rules are: the source read segment number remains unchanged, the source read segment starting position is set to a position that is (L+1) (backward) away from the original source read segment starting position, and the source read segment ending position remains unchanged.

[0115] c. The path unit objects generated by the above clustering and mapping may include path unit objects with empty path unit sequences. Therefore, it is necessary to perform step 5) to merge equivalent tangent point IDs on the obtained objects. The path unit objects obtained after merging equivalent tangent point IDs are the path unit objects after left tangent alignment finally generated in step 6).

[0116] 7) Perform right tangent point alignment on all path unit objects generated in the previous step. From steps 1)-2), it can be obtained that the path bifurcation point in the final variation path graph must correspond to a tangent point. Therefore, the first character on the left side of a tangent point determines the unique path between two adjacent tangent points in the variation path graph. Based on this principle, the path unit can be right-aligned. The right alignment of the path unit can split the fused path unit. The existence of the fused path unit is due to the truncation effect on the left and right sides of the read segment, which causes some tangent points within the length (K-1) on the left and right sides of the read segment to be unable to be marked by steps 1)-2), so that the path unit obtained in steps 3)-4) is not the smallest path unit, but a fusion of several path units that cross the unmarked tangent points. The steps for right-aligning the path unit are as follows:

[0117] a. For all path unit objects, cluster them according to their (end cut point ID, the last character of the path unit sequence) tuples;

[0118] b. Map each cluster of path unit objects to a new cluster of path unit objects as follows:

[0119] 1. Sort the path unit objects in the cluster according to the length of the path unit sequence. The path unit object with the smallest length is A, and its length is L.

[0120] 2. Map A to a new path unit object. The information of the new object is consistent with A, and only the source read segment position table of the object is changed and updated. The source read segment position table of the new object contains:

[0121] i) All records in the source read position table of A;

[0122] ii) Records mapped from records in the source read position table of other path units in the cluster. The mapping rule is: the source read segment number remains unchanged, the source read segment start position is set to a position L (forward) away from the source read segment end position, and the source read segment end position remains unchanged.

[0123] 3. Map the non-A path unit object in the cluster to a new path unit object. The starting cut point ID of the new object remains unchanged. The ending cut point ID of the new object is set to the starting cut point ID of A. The path unit sequence of the new object is set to the sequence remaining after deleting the path unit sequence of A from the right side of the path unit sequence of the original object (the sequence may be empty). The records of the source read position table of the new object are mapped from the records of the source read position table of the original object. The mapping rules are: the source read segment number remains unchanged, the source read segment starting position remains unchanged, and the source read segment ending position is set to a position that is (L+1) (forward) away from the original source read segment ending position.

[0124] c. The path unit objects generated by the above clustering and mapping may contain path unit objects with empty path unit sequences, so it is necessary to perform the operation of merging equivalent tangent point IDs in step 5) on the obtained objects. The path unit objects obtained after merging equivalent tangent point IDs are the path unit objects after right tangent point alignment finally generated in step 7).

[0125] 8) Repeat steps 6)-7) until each cluster contains only one path unit object during the execution of steps 6) and 7).

[0126] 9) Based on the path unit object generated in step 8), generate a "path segment internal tangent point set", where the path segment internal tangent point points are tangent points that are connected to only one path unit (i.e., not bifurcated) both forward and backward. The steps for generating the path segment internal tangent point set are as follows:

[0127] a. Cluster all path unit objects according to the starting cut point ID, and record the starting cut point ID of the cluster containing only one path unit object, which is recorded as set A.

[0128] b. Cluster all path unit objects according to the end cut point ID, and record the end cut point ID of the cluster containing only one path unit object, which is recorded as set B.

[0129] c. Take the intersection of sets A and B to obtain the internal tangent point set of the path segment.

[0130] 10) Map the path unit object generated in step 8) to an unconnected "path segment" object. Each path segment object contains the following information:

[0131] a. Starting point ID (equal to the starting point ID of the path unit object);

[0132] b. Whether the starting point of contact is an internal point of contact of the path segment (if the starting point of contact ID is in the internal point of contact set of the path segment generated in step 9), then yes, otherwise no);

[0133] c. End cut point ID (equal to the end cut point ID of the path unit object);

[0134] d. Whether the end cut point is an internal cut point of the path segment (if the end cut point ID is in the internal cut point set of the path segment generated in step 9), then yes, otherwise no);

[0135] e. Path segment length (set to the length of the path unit sequence);

[0136] f. Source read mapping table. A record in each path unit source read position table is mapped to a record in the path fragment source read mapping table, which contains the following information:

[0137] The path segment starting position (set to 1);

[0138] The path segment end position (set to the length of the path unit sequence);

[0139] Read segment ID (equal to the source read segment ID recorded in the path unit source read segment position table);

[0140] The starting position of the read segment (equal to the source read segment starting position recorded in the path unit source read segment position table);

[0141] The end position of the read segment (equal to the source read segment end position recorded in the path unit source read segment position table).

[0142] 11) Connect the unconnected path segment objects generated in step 10), as follows:

[0143] a. Find a pair of connectable path segment objects, denoted as A and B, that meet the following conditions: the end cut point ID of A is equal to the start cut point ID of B.

[0144] b. Combine A and B into a new path segment object C. After the combination, A and B no longer exist. C can continue to connect with other path segment objects if the conditions are met. The combination rules are as follows:

[0145] The starting cut point ID of C is set to the starting cut point ID of A;

[0146] C's "Is the starting tangent point an internal tangent point of the path segment" is set to A's "Is the starting tangent point an internal tangent point of the path segment"

[0147] "Internal tangent point of segment";

[0148] The end point ID of C is set to the end point ID of B;

[0149] C's "Is the end point of the path segment internal?" is set to B's "Is the end point of the path segment internal?"

[0150] "Internal tangent point of segment";

[0151] The path segment length of C is set to the sum of the path segment length of A and the path segment length of B;

[0152] The source read mapping table of C is generated by merging the source read mapping tables of A and B. The steps are as follows:

[0153] i) Since the B path segment is connected after the A path segment, first increase the “read segment start position” and “read segment end position” of each record in the source read segment mapping table of the B path segment by the “path segment length” of the A path segment;

[0154] ii) Take the source read mapping tables of the two path segments A and B, and obtain the union of their records. In this set, if two records (referred to as U and V) meet the conditions "U comes from the source read mapping table of A; V comes from the source read mapping table of B; the read IDs of U and V are the same; the end position of the read of U is equal to the start position of the read of V minus 1", then combine these two records into a new record (referred to as R), so that "the path segment start position of R is equal to the path segment start position of U; the path segment end position of R is equal to the path segment end position of V; the read ID of R is equal to the common read ID of U and V; the read start position of R is equal to the read start position of U; the read end position of R is equal to the read end position of V";

[0155] iii) The record set generated by the above steps constitutes the source read segment mapping table of C.

[0156] c. Repeat the above steps until no connectable path segment objects can be found. The resulting path segment object set is the fully connected path segment object set.

[0157] 12) The fully connected set of path segment objects generated by the above steps (each path segment object contains the cut point ID information on its left and right sides, as well as which part of the path segment sequence is covered by which fragment of the original sequencing data read sequence) constitutes a variation path graph representation of the original sequencing data.

[0158] In view of the massive amount of raw data, the present invention provides another specific method for generating a variation path map of sequencing data. This method divides the massive amount of raw sequencing data into different parts, constructs a variation path map for each part of the raw sequencing data, and then merges the obtained variation path maps to obtain a variation path map of all sequencing data or specific sequencing data, comprising the following steps:

[0159] 1) Splitting the massive amount of raw sequencing data into different parts according to samples or groups, and generating a variation path map of each part of the raw sequencing data according to the method described in claim 2 above;

[0160] 2) Merging the obtained variation path graphs in pairs to obtain the variation path graph of all original sequencing data; or merging the variation path graphs of specific samples or groups to obtain the variation path graph of a part of specific sequencing data, which is used to compare the differences between specific samples or specific groups. The steps of merging two variation path graphs are as follows:

[0161] a. Give each position of the path segment in the two variant path graphs (denoted as G and H) a unique position number, so that:

[0162] In the same variation path graph, all positions of all path segments have different numbers.

[0163] In the same variation path diagram, the numbers from the leftmost to the rightmost position of the same path segment increase continuously.

[0164] In the same variation path graph, the numbers of any two positions belonging to two different path segments are not adjacent (the difference is not equal to 1).

[0165] b. Treat the mutation path graphs G and H as a read segment respectively. G has a read segment ID "1" and H has a read segment ID "2". Scan all strings of length (K-1) in G and H to generate (K-1) string objects. These strings may start from any position in any path segment and cross the boundaries of the path segment in any path. Each (K-1) string object contains the following information:

[0166] The scanned sequence string is of length (K-1).

[0167] The characters on the left side of the string (if the string comes from the leftmost side of a path segment that is not connected to other path segments on the left, the characters on its left side are recorded as '|'; if the string comes from the leftmost side of a path segment that may be connected to multiple other path segments on the left side, multiple (K-1) string objects are generated, and each (K-1) string object records one possibility).

[0168] The characters on the right side of the string (if the string comes from the rightmost side of a path segment that is not connected to other path segments on the right side, the characters on the right side are recorded as '|'; if the string comes from the rightmost side of a path segment that may be connected to multiple other path segments on the right side, multiple (K-1) string objects are generated, and each (K-1) string object records one possibility).

[0169] The source read segment number of the (K-1) string ("1" or "2").

[0170] The starting position number of the (K-1) string in the source read segment.

[0171] The ending position number of the (K-1) string in the source read segment.

[0172] c. According to the method similar to the aforementioned claim 2, based on the above (K-1) string objects, a synthetic variation path graph can be generated, and the main changes in the method are:

[0173] The (K-1) string object no longer records the full sequence information of its source read segment, but records the end position of the (K-1) string on the read segment. This change is because the (K-1) string may cross the boundary of the path segment, and the position numbers in different path segments are not continuous, so its end position cannot be calculated based on the starting position of the K-1 string.

[0174] When generating the "cut point-read segment position" object, the end position information of the (K-1) string recorded by the (K-1) string object on the source read segment is directly used, rather than calculating the end position based on its starting position. At the same time, the cut point-read segment position object needs to record the characters on both sides of the cut point (which can be obtained based on the information recorded by the (K-1) string object), but does not record the full sequence of the source read segment.

[0175] When generating a 'pathunit-read position' object, the path unit sequence is not recorded, but the first and last characters of the path unit are recorded.

[0176] Likewise, when a PathUnit object is generated, the path unit sequence is not recorded, but the first and last characters of the path unit are recorded.

[0177] Since there are only two reads, there will be no "fused path unit", so all steps related to "left cut point alignment" and "right cut point alignment" are skipped.

[0178] The path units generated by the above method must correspond to some split of the path units in the original variant path graphs G and H. Therefore, when generating unconnected "path segments", the length of the "path unit sequence" of the path unit can be calculated based on the distance between the "source read starting position" and the "source read ending position" of any record in its source read position table (at this time, these two positions must be located in the same path segment of G or H, so the difference in position numbers is equal to the sequence length).

[0179] The synthetic variation path graph (a set of fully connected path segment objects) generated by the above method must have the starting and ending positions of the path segments corresponding to the starting or ending positions of a path segment in the original variation path graph G or H, that is, the path segment of the synthetic variation path graph must be the split of the path segments in G and H. Based on this correspondence, the source read segment mapping table of the path segments in G and H can be converted to the source read segment mapping table of the path segments in the synthetic variation path graph. The conversion method is an intuitive segment splitting method, which will not be repeated here.

[0180] d. According to the above method, the two variation path graphs G and H can be merged into a synthetic variation path graph.

[0181] The above is not a limitation on the specific implementation methods of this patent. It should be pointed out that for ordinary technicians in this technical field, several changes, modifications, additions or substitutions can be made without departing from the essential scope of the present invention, and these improvements and modifications should also be regarded as the protection scope of the present invention.

Claims

1. A method for generating a variation path graph of sequencing data, characterized in that: By identifying and tracking the cut points, the association between "the bifurcation point of the pathway segment or pathway unit in the variation pathway map", "the position inside the pathway segment sequence in the variation pathway map corresponding to the starting point or end point of the original sequencing read, i.e., the starting point or end point of the pathway unit" and "the position coordinates of the sequence in the original sequencing read" is established; The steps of obtaining a variation path map include:

1. identifying the position of the cut point in the original read sequence; 2. splitting the original read sequence into small fragments marked by the cut points at the beginning and end; 3. grouping these small fragments according to the first and last cut points and the sequences in between to generate path units; 4. connecting the non-branched path units head to tail to form path fragments and obtain a variation path map; the variation path map includes path fragments and cut points; the path fragments include path units; the beginning and end of the path fragments and the path units are marked by cut points; if the end cut point of a path fragment is the same as the beginning cut point of another path fragment, they are connected from left to right; the beginning and end of a path fragment can be connected to multiple path fragments or not connected to the path fragment; it is called path bifurcation; the path fragment is connected by path units; the beginning and end of the path unit are also marked by cut points, and the end cut point of one path unit is the same as the beginning cut point of another path unit The first cut points of the two reads are the same, then they are connected from left to right; the path units inside the path fragment are connected one-to-one in order without bifurcation; the first and last path units of the path fragment are either not connected to other path units, or connected to multiple path units belonging to other path fragments respectively; a read segment in the original sequencing data covers a unique path in the variation path graph; the starting and ending positions of each read segment must correspond to a cut point; if there is an identical subsequence of length "K-1" in two read segments, but the left or right characters adjacent to the identical subsequence are different, then the left or right position of the identical subsequence of the two read segments corresponds to the same cut point in the variation path graph, that is, the cut point at the bifurcation of the path; if there is an identical subsequence of length "K-1" in two read segments, then the identical subsequence covers the same and unique path in the variation path graph.

2. The method for generating a variation path graph of sequencing data according to claim 1, characterized in that: The specific steps are as follows: 1) First, scan the read sequences of all the original sequencing data, generate all possible subsequences of length "K-1", and record the characters on the left and right sides of the subsequences. Then, cluster these subsequences of length "K-1", and identify the position of the cut point in the original read sequence according to the possibility of the characters on the left or right side of each type of subsequence of length "K-1"; 2) Cluster all the "K-1" string objects generated in step 1) according to their "K-1" length strings, assign a unique cluster ID to each cluster of "K-1" string objects, where the ID is a positive integer, and then analyze the "K-1" string objects of each cluster to obtain the cut point-read segment position object; 3) clustering all the cut point-read position objects generated in step 2) according to their source read numbers; Then the cut point-read position object of each cluster is analyzed to generate a pathway unit-read position object; 4) clustering all the path unit-read position objects generated in step 3) according to the triplets of the start cut point ID, the end cut point ID, and the path unit sequence; and then mapping each cluster of path unit-read position objects into a path unit object; 5) For all the path unit objects generated in the previous step, perform the operation of merging equivalent tangent point IDs; 6) Perform left tangent point alignment operation and right tangent point alignment operation; 7) Repeat step 6) until during the execution of step 6), each cluster contains only one path unit object; 8) Based on the path unit object generated in step 7), a set of internal tangent points of the path segment is generated, wherein the internal tangent points of the path segment are only connected to one path unit forward and backward, i.e., non-branching tangent points; 9) Mapping the path unit object generated in step 7) into an unconnected path segment object; 10) connecting the unconnected path segment objects generated in step 9); 11) The fully connected set of path fragment objects generated by the above steps constitutes a variation path graph representation of the original sequencing data; the path fragment object includes the cut point ID information on the left and right sides, and which part of the path fragment sequence is covered by which fragment of the original sequencing data read sequence.

3. A method for generating a variation path graph of sequencing data, characterized in that: The original sequencing data is divided into different parts, and a variation path map is constructed for each part of the original sequencing data according to the method of claim 1, and then the obtained variation path maps are merged to obtain a variation path map of all sequencing data or specific sequencing data.

4. The method for generating a variation path graph of sequencing data according to claim 3, characterized in that: The specific steps of merging the two variation path graphs are as follows: 1) Give each position of the path segment in the two variant path graphs a unique position number; 2) Treat the two variation path graphs as a read segment respectively, scan all the strings of length "K-1" in the two read segments, and generate a "K-1" string object. The string may start from any position in any path segment and cross the boundary of the path segment with any path; the "K-1" string object no longer records the full sequence information of its source read segment, but records the end position of the "K-1" string on the read segment; 3) Clustering the "K-1" string objects generated in step 2) according to their "K-1" length character strings, assigning a unique cluster ID to each cluster's "K-1" string object, where the ID is a positive integer, and then analyzing each cluster's "K-1" string object respectively, directly using the end position information of the "K-1" string on the source read segment recorded by the "K-1" string object, rather than calculating the end position based on its starting position, to obtain a cut-point-read segment position object; at the same time, the cut-point-read segment position object needs to record the characters on both sides of the cut-point, but does not record the full sequence of the source read segment, and the characters on both sides of the cut-point can be obtained based on the information recorded by the "K-1" string object; 4) clustering all the cut point-read position objects generated in step 3) according to their source read numbers; then analyzing the cut point-read position objects of each cluster respectively, performing operations to generate path unit-read position objects; when generating the path unit-read position objects, the path unit sequence is not recorded, but the first and last characters of the path unit are recorded; 5) clustering all the path unit-read position objects generated in step 4) according to the four-tuple of the starting cut point ID, the ending cut point ID, and the first and last characters of the path unit; then mapping each cluster of path unit-read position objects into a path unit object; when generating the path unit object, the path unit sequence is not recorded, but the first and last characters of the path unit are recorded; 6) For all the path unit objects generated in the previous step, perform the operation of merging equivalent tangent point IDs; 7) Based on the path unit object generated in step 6), a set of internal tangent points of the path segment is generated, wherein the internal tangent points of the path segment are tangent points that are connected to only one path unit both forward and backward, i.e., tangent points that are not bifurcated; 8) Mapping the path unit object generated in step 6) to an unconnected path segment object, where the length of the path unit sequence of the path unit can be calculated based on the distance between the source read start position and the source read end position of any record in its source read segment position table; 9) connecting the unconnected path segment objects generated in step 8); 10) The start and end positions of the fully connected path segment object set generated by the above steps must correspond to the start position or end position of a path segment in the original variation path graph, that is, the path segment of the synthesized variation path graph must be the split of the path segment in the original variation path graph. Based on this correspondence, the source read segment mapping table of the path segment in the original variation path graph can be converted into the source read segment mapping table of the path segment in the synthesized variation path graph; 11) The set of path segment objects generated by the above steps is the result of merging the two variation path graphs, forming a variation path graph representation of the merged original sequencing data; the path segment object includes the cut point ID information on its left and right sides, and which part of the path segment sequence is covered by which fragment of the merged original sequencing data read sequence.

Images