A method for subgenomic identification in a polyploid genome assembly process

Abstract

The application discloses a method for identifying sub-genome in a polyploid genome assembly process, and optimizes results by connecting any Query and Subject HSP in series through a modified dynamic programming algorithm, wherein the core innovation point is that the coordinates of the HSP of the Subject are replaced by the coordinates mapped on the Query, so that the influence of the rearrangement event of the Subject on the HSP chain analysis result can be maximally excluded, and the application relates to the technical field of genome assembly. In the method for identifying the sub-genome in the polyploid genome assembly process, only the 1vs 1 alignment result is used in the analysis process in the self-alignment elimination method, and logical loopholes (the 'bridge' event cannot be recognized) exist. Since the modified dynamic programming algorithm used by the application can ignore the influence of the rearrangement in the Subject on the HSP chain, the application can use the modified any two Subjects similar to the Query to do any connection to form a new Subject, and the possible 'bridge' event can be recognized.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of genome assembly, in particular to a method for identifying subgenomes in the process of assembling polyploid genomes. BACKGROUND

[0002] Polyploid genomes refer to the presence of three or more complete sets of chromosomes in the cells of an organism. This is different from the usual diploid organisms (such as humans), whose cells contain only two sets of chromosomes. Polyploidy is particularly common in plants, but also exists in some animals.

[0003] Main types of polyploid genomes:

[0004] 1. Homoploids: multiple sets of chromosomes from the same species.

[0005] 2. Allopolyploids: multiple sets of chromosomes from different but related species.

[0006] Characteristics of sequence similarity between polyploid chromosomes:

[0007] 1. High similarity: Chromosome sets in homoploids usually have high sequence similarity. This similarity can lead to increased complexity in genome assembly and analysis.

[0008] 2. Homologous regions: There are a large number of homologous regions between different chromosome sets. These regions may be almost completely identical in sequence, especially in recently formed polyploids.

[0009] 3. Gradual differentiation: Different chromosome sets may gradually differentiate over time. This differentiation can lead to a decrease in sequence similarity and the generation of subgenome-specific variations.

[0010] 4. Gene redundancy: Multiple almost identical gene copies exist in different chromosome sets. This redundancy can lead to complex regulation of gene expression and function.

[0011] 5. Structural variations: Rearrangement, insertion, deletion, etc. structural variations may occur between chromosomes. These variations can cause differences in local sequence similarity.

[0012] 6. Homologous recombination: Homologous recombination may occur between highly similar chromosomes. This can lead to gene conversion and sequence assimilation.

[0013] 7. Subgenome preference: In some polyploids, there may be a phenomenon of preferential preservation or expression of a certain subgenome. This can lead to different degrees of sequence conservation between subgenomes.

[0014] 8. Repeat distribution: The distribution of transposable elements and other repeats can be different among different chromosomes. This can cause differences in local sequence similarity.

[0015] 9. Functional redundancy and new functions: Highly similar sequences can lead to functional redundancy of genes. At the same time, it also provides opportunities for the differentiation of gene functions and the generation of new functions.

[0016] 10. Sequence loss: After polyploidization, loss of genes or chromosomal fragments can occur. This phenomenon can cause originally similar regions to become different.

[0017] In summary, there is a certain similarity between the sequences of polyploid chromosome sister chromatids, but due to independent evolution and other conditions, there are large differences in gene rearrangement, insertion and deletion, etc. Understanding these characteristics is of great significance for accurately analyzing polyploid genome structure, function and evolution, but also makes the assembly and analysis of polyploid genomes complex.

[0018] After genome assembly analysis, a consistent sequence representing a haploid genotype is expected. After the preliminary assembly of the polyploid genome data of a species by assembly software, the contig or scaffold version of the genome is obtained, which can be considered as a fragmented sequence of the genome chromosome. These contigs must contain homologous sequences of different homologous groups or different subgenomes. This will result in a large number of homologous sequences representing the same chromosomal region in these homologous groups, which have certain similarities and certain differences, deletions and rearrangements between each other. The existence of these redundant homologous sequences will affect subsequent analysis (such as an increase in the number of genes, incorrect chromosome localization), and need to be removed as much as possible.

[0019] The conventional methods are:

[0020] 1. Reference sequence alignment method:

[0021] The reference sequence alignment method generally uses a close relative species reference sequence. These contigs are aligned to the reference sequence using high-mismatch-tolerant alignment tools such as Minimap2 and Lastz. Then a threshold is set. Sequences below the threshold will be excluded, and the remaining sequences will be retained as the result.

[0022] 2. Self-alignment elimination method:

[0023] All contigs are paired using an alignment tool. Then, the local alignment results are concatenated using the Lastz Chain algorithm. The longer of the two best-aligned sequences is retained, and the shorter one is removed as a subgenome. Of course, some tools use coverage as the criterion and remove homologous contigs with low coverage as subgenomes, such as HaploMerger2, Redundans, and purge_haplogs.

[0024] 3. Spectrum Comparison Method:

[0025] Using BioNano or HiC sequencing results, all contigs are compared with the map. Within the same comparison region, the contig with the highest score is retained, and the other contigs in the same region are removed as subgenomes.

[0026] However, the alignment tools used in the first two methods, such as LastZ and Minimap2, employ the Chain algorithm primarily designed to obtain collinearity results between genomes. For reverse complementarity or skip rearrangements, these algorithms, by default, do not recognize the alignment results, affecting the final removal of subgenomes (e.g., Figure 1 ).

[0027] The map comparison method also has certain limitations:

[0028] 1) The sister chromatid used to construct the map may not be the same chromatid obtained from the actual assembly.

[0029] 2) For small individual samples, all samples used to construct the atlas may not be the same samples used for sequencing, and sister chromatids may undergo homologous recombination during sexual reproduction;

[0030] 3) The assembled Contig is a consistent sequence, meaning that its assembly result may be a mixture of multiple homologous chromosome types, and there may be misspellings within it.

[0031] Therefore, during polyploid assembly, the alignment of the spectra and contigs inevitably results in numerous rearrangements and differences. The conventional Chain algorithm's neglect of large-scale rearrangements leads to a large number of subsequent analysis errors, requiring manual correction.

[0032] For these challenges, the common approach is to use dynamic programming algorithm and its derived heuristic or greedy optimization algorithm to chain the HSPs (High-scoring Segment Pair) from local alignments (which can be Contig-to-Contig alignment, Contig-to-reference alignment, or Contig-to-genetic map alignment) and calculate their similarity (the length of the chained alignment / the total length of the Contig) to determine and remove the redundant homologous Contig sequences. Examples Figure 2 are shown.

[0033] a) Suppose there are two sequences Refl (5400bp) and Ref2 (7600bp), Refl and Ref2 are highly homologous, and Refl can be completely contained in Ref2 without considering rearrangement. Figure 2 The regions with the same color represent the same sequence, the arrow from left to right direction represents 5'→3', and the arrow from right to left direction represents 3'←5'. The box represents the sequence unique to Ref2, which has no homology with Refl. The gray region is inverted, and the black region is inverted and rearranged.

[0034] b) The vertical coordinate represents the alignment coordinate of Refl, and the horizontal coordinate represents the alignment coordinate of Ref2. The HSPs obtained by aligning Refl and Ref2 are sorted according to the coordinates to construct a matrix. If two regions can be aligned, the alignment score is placed in the corresponding position of the matrix.

[0035] c) Use dynamic programming to chain the HSPs.

[0036] Principle:

[0037] 1. Use dynamic programming algorithm to connect similar HSPs into longer sequences.

[0038] 2. Consider the distance, direction, and score between HSPs to find the highest scoring position.

[0039] Algorithm steps:

[0040] 1. Sort all HSPs according to their starting positions on the reference sequence.

[0041] 2. Traverse the sorted HSPs from left to right: a. For each HSP, try to connect it with the previous HSPs (left or top-left). b. Use dynamic programming to calculate the best connection method.

[0042] 3. Dynamic programming score calculation: score(i) = max{score(j) + match_score(i) - gap_penalty(i, j)} where:

[0043] i is the current HSP, j is the previous HSP, match_score(i) is the score of the current HSP, gap_penalty(i,j) is the gap penalty when connecting i and j (for simplicity, all gap_penalty in this example is 0);

[0044] 4. Find the matrix position with the highest score

[0045] 4. Backtrack the score matrix to find the best chain path. (Red and green arrow positions)

[0046] 5. According to the set threshold, filter out the chains with lower scores.

[0047] 6. Output the final chain result.

[0048] 7. Figure 2 The HSP represented by the medium gray arrow has part of the Ref1 coordinates (4500-5400) that are repeatedly calculated, so only one with the highest score is retained in the two regions for chain calculation.

[0049] The core idea of this algorithm is to maximize the overall alignment score while maintaining sequence colinearity.

[0050] It can be found that under the premise of maintaining sequence colinearity, Ref1 can be aligned with Ref2 for 5400bp using conventional algorithms, but the optimal calculation result is only 4300bp. Assuming that the threshold for merging homologous sequences is 85% similarity, Ref1 sequence will be missed, resulting in inaccurate calculation results.

[0051] In addition, for the self-alignment elimination method, since the current algorithm focuses on 1vs1 similarity alignment and processing between sequences, but for Figure 3 As shown in the case (the Contig in the dashed box is completely contained by the other two Contigs, but no Contig can meet the screening threshold in 1vs1 alignment, which we call "bridging" events), there is currently no algorithm that can automatically handle this logical loophole, relying on manual inspection, which is time-consuming and laborious.

[0052] To solve this problem, the following technical solutions are adopted:

[0053] 1. Manual inspection of alignment results, manually checking the alignment results of each allele, and manually deleting homologous sequences (reference comparison file CN111445948A CN111584004A);

[0054] 2. Local modification method, for the alignment results of genetic map and reference sequence. If multiple Contigs can be aligned to the same reference sequence or a region of the genetic map, the sequencing Reads constituting these Contigs are collected, then local reassembly or clustering is performed, and then the reassembled Contigs are aligned to the map again until a Contig satisfying the threshold is obtained, and the other Contigs are discarded (for comparison file CN109326323B CN110020726A CN112289382A);

[0055] 3. Sample pretreatment method, construct SNP linkage groups before sequencing. According to the linkage group, select the original data and construct the genetic map, and only analyze the haploid data data satisfying the same SNP linkage group, so as to solve the problem of polyploidy from the source (for comparison file CN112289382B, CN111816248A).

[0056] The above three methods are time-consuming and laborious, and difficult to repeat. At the same time, it is inevitable that there is a certain difference between the genetic map and the Contig, and local modification may not be able to obtain the optimal result, and excessive modification may damage the gene structure, resulting in a series of subsequent problems. SUMMARY

[0057] (1) Technical problems solved

[0058] In view of the deficiencies of the prior art, the present application provides a new subgenome identification method in the process of assembling polyploid genome based on a modified dynamic programming algorithm, which can eliminate the influence of rearrangement events on the analysis results. Compared with the existing method, this method can more accurately identify subgenomes in polyploidy, and can perfectly solve the logical loophole of "bridging" events.

[0059] (2) Technical solutions

[0060] In order to achieve the above purpose, the present application is implemented by the following technical solutions: a subgenome identification method in the process of assembling polyploid genome, specifically comprising the following steps:

[0061] S1, obtaining the alignment result of any Contig (Query), for any Query, obtaining the alignment coordinate matrix of any target Subject;

[0062] S2, the matrix is characterized in that one-dimensional coordinates are the alignment coordinates of Contig, and the other dimension coordinate system is the coordinates of Subject mapped to Query;

[0063] S3, for any one cell G(i, j), score(i, j) = max{score(m, n) + match_score(i, j) - gap_penalty(m, n)(i, j)}, where i > 0, j > 0, m <= i, n <= j;

[0064] S4, scan the whole matrix to find the highest score, and use dynamic programming algorithm to calculate the sum of base coordinates involved in the path of the highest score X;

[0065] S5, calculate the ratio of X and the length of Query sequence, if it is greater than a set threshold, then mark Query as the homologous sequence of Subject and remove it;

[0066] S6, for any one Contig set S, the above steps S1-S5 algorithm can split Contig into two sets, S 去除 and S 保留 ;

[0067] S7, to eliminate the logical loophole appeared in self-comparison, for any one C in S 保留 , query the previous dynamic programming results to find the first two Contigs C x and C y contained in S 保留 with the highest coordinate sum of C, concatenate the sequences of C x and C y to get the sequence C connect , compare C connect with C, then perform steps S1-S6 on them and update S 保留 ;

[0068] S8, traverse S 保留 until the end.

[0069] Preferably, assuming there are two sequences Ref1 and Ref2, Ref1 and Ref2 are highly homologous, and Ref1 can be completely contained in Ref2 without considering rearrangement.

[0070] Preferably, the number of base pairs of Ref1 is 5400 bp, and the number of base pairs of Ref2 is 7600 bp.

[0071] Preferably, the vertical coordinate is the Ref1 alignment coordinate, the horizontal coordinate represents the Ref2 HSP mapping to the Ref1 alignment coordinate, and the matrix is constructed according to the new coordinate sorting, if two regions can be aligned, the alignment score is placed in the corresponding position of the matrix.

[0072] Preferably, the HSPs are concatenated using dynamic programming, and the formula is as follows:

[0073] score(i) = max{score(j) + match_score(i) - gap_penalty(i,j)};

[0074] Where the coordinate of j cell is less than or equal to the coordinate of i cell.

[0075] Where: i is the current HSP, j is the previous HSP, match_score(i) is the score of the current HSP, and gap_penalty(i,j) is the gap penalty when connecting i and j.

[0076] Preferably, the result shows that the region of the alignment obtained by the Chain is 5200bp, which is significantly better than the 4300bp obtained by the conventional algorithm.

[0077] (Three) beneficial effects

[0078] The application provides a method for identifying subgenomes in a polyploid genome assembly process.

[0079] (1) The method for identifying subgenomes in a polyploid genome assembly process concatenates any Query and Subject HSP to obtain an optimal result through a modified dynamic programming algorithm, and the core innovation point is that the coordinates of the HSP of the Subject are replaced with coordinates mapped to the Query, so that the influence of the rearrangement event of the Subject on the analysis result of the HSP concatenation (Chain) can be maximally excluded.

[0080] (2) The method for identifying subgenomes in a polyploid genome assembly process can concatenate HSPs of any source through the use of a modified programming algorithm, which can represent sequence alignment coordinates, genetic map or SSR molecular marker map alignment coordinates, enzyme digestion map or HiC alignment matrix, etc. As long as a tensor or vector related to similarity scoring and sorting is involved, the method can be applied. Therefore, the protection scope needs to be expanded. The examples given only use the start coordinates of the sequence alignment result as an example for easy understanding, but in fact, the method can be applied to any vector or tensor representing a similarity alignment system.

[0081] (3) The method accepts a sequence formed by randomly concatenating multiple Query sequences as input, and since the method can exclude the influence of rearrangement events on the analysis result, it can automatically and perfectly solve the "bridging" problem. BRIEF DESCRIPTION OF DRAWINGS

[0082] Figure 1 Figure for the shortcomings of the existing LastZ and Minimap2 alignment result processing;

[0083] Figure 2 Figure for the dynamic programming Chain algorithm under the conventional framework;

[0084] Figure 3 Figure for the logical loopholes in the current homologous Contig screening process;

[0085] Figure 4 Figure for the HSP after the series connection of different methods of the present application;

[0086] Figure 5 Figure for the framework of the Chain algorithm of the present application;

[0087] Figure 6 Figure for the first demonstration of the Chain algorithm of the present application;

[0088] Figure 7 Figure for the second demonstration of the Chain algorithm of the present application;

[0089] Figure 8 Figure for the first demonstration of the Chain algorithm of the present application;

[0090] Figure 9 Figure for the second demonstration of the Chain algorithm of the present application. DETAILED DESCRIPTION

[0091] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0092] Please refer to Figures 1-9 The embodiments of the present application provide a technical solution: a method for identifying subgenomes in the assembly of polyploid genomes, which specifically includes the following embodiments:

[0093] A method for identifying subgenomes in the assembly of polyploid genomes, which specifically includes the following steps:

[0094] S1, obtaining the alignment result of any Contig (Query), and obtaining the alignment coordinate matrix of any Query and any target Subject;

[0095] S2, the matrix is characterized in that its one-dimensional coordinate is the alignment coordinate of Contig, and the coordinate system of the other dimension is the coordinate of Subject mapped onto Query;

[0096] S3, for any one cell i,j, score(i,j)=max{score(m,n)+match_score(I,j)-gap_penalty({m,n}(i,j))}, wherein i>=0,j>=0,m<=i,n<=j;

[0097] S4, scanning the whole matrix to find the highest score value, and using the dynamic programming algorithm to calculate the base coordinate sum X involved in the path of the highest value;

[0098] S5, calculating the ratio of X and Query, if greater than the set threshold, then mark Query as the homologous sequence of Subject, and remove it;

[0099] S6, for any one Contig set S, the above steps S1-S5 algorithm can split Contig into two sets, S 去除 and S 保留 ;

[0100] S7, for the logical loophole appeared in the self-alignment elimination method, from S 保留 , any one C is operated, the previous dynamic programming result is inquired, the first two Contigs C 保留 contained by S x and C y are found, C x and C y are concatenated to obtain sequence C connect , C connect is aligned with C, then steps S1-S6 are performed on C and C, and S 保留 is updated;

[0101] S8, S 保留 is traversed until the end.

[0102] In the embodiment of the application, it is assumed that there are two sequences Ref1 and Ref2, Ref1 and Ref2 are highly homologous, Ref1 can be completely contained by Ref2 without considering rearrangement, Figure 4The regions with the same color represent the same sequence, the arrow from left to right direction represents 5'→3', from right to left represents 3'←5' direction, the block represents the sequence unique to Ref2, and Ref1 has no homology. The light gray region is flipped, the gray region is jump-flipped rearrangement, the number of Ref1 base pairs is 5400bp, and the number of Ref2 base pairs is 7600bp, the vertical coordinate is the standard Ref1 alignment coordinate, the horizontal coordinate represents the alignment coordinate of Ref2 HSP mapped to Ref1, and the matrix is constructed according to the new coordinate sorting. If two regions can be aligned, the alignment score is placed in the corresponding position of the matrix, the HSP is connected by using dynamic programming, and the calculation formula is as follows:

[0103] score(i)=max{score(j)+match_score(i)-gap_penalty(i,j)};

[0104] Wherein: i is the current HSP, j is the previous HSP, match_score(i) is the score of the current HSP, and gap_penalty(i,j) is the gap penalty when connecting i and j (for simplicity, all gap_penalty in this example is 0).

[0105] It can be found that the region aligned by this Chain can reach 5200bp, which is significantly better than 4300bp of the conventional algorithm.

[0106] The comparison results of the sample data of the polyploid chelicerate using the method and the conventional method are shown in Table 1.

[0107] Table 1 statistical table of examples

[0108]

[0109]

[0110] As Figure 1A is the vulnerability of the algorithm related to the Lastz algorithm. Assuming that the color box is the result of the reverse complementary alignment, Lastz ignores the reverse complementary alignment result that appears intermittently in the process of concatenating HSPs (High-scoring Segment Pair). B is the common vulnerability of Minimap and Lastz. For the jump rearrangement that occurs universally in sub-chromosomes, the algorithm of this software will ignore it by default. In addition, these tools need to be annotated with repeat sequences in advance to prevent repeat sequences from affecting the analysis results. Finally, for large-scale rearrangements, such as the deletion of a part of the whole chromosome or jumping to other positions, these algorithms will consider them to have lost homology, and the final alignment result will not be recognized by the Chain algorithm, thereby affecting the analysis accuracy.

[0111] Figure 5 The left graph is the LastZ Chain algorithm, and the right graph is the method of this application (the result of the LastZ algorithm is retained for easy checking and comparison).

[0112] Figure 6 The middle two sequences are the alignment results of the two sequences. Assuming that there are Ref1 (5400bp) and Ref2 (7600bp) two sequences, Ref1 and Ref2 are highly homologous, and Ref1 can be completely contained in Ref2 without considering rearrangement. The same color area in the figure represents the same sequence, the arrow from left to right direction represents 5'→3', and the arrow from right to left direction represents 3'←5'. The box represents the sequence unique to Ref2, which has no homology with Ref1. The gray-white area is inverted, and the gray area is jump-inverted rearrangement.

[0113] Meanwhile, the contents not described in detail in the specification are all prior art known to those skilled in the art.

[0114] It should be noted that, in this text, relational terms such as first and second are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "contain" or any other variant thereof are intended to cover non-exclusive inclusion, so that the process, method, article or equipment including a series of elements not only includes those elements, but also includes other elements not explicitly listed or inherent to such process, method, article or equipment.

[0115] Although embodiments of the present application have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made thereto without departing from the principles and spirit of the present application, and the scope of the present application is defined by the appended claims and their equivalents.

Claims

1. A method of subgenomic identification in a polyploid genome assembly process, characterized by: Specifically comprising the following steps: S1, obtaining the alignment result of any Contig, for any Query, obtaining the alignment coordinate matrix of any target Subject; S2, the matrix is characterized in that one-dimensional coordinate is the alignment coordinate of Contig, the other dimension coordinate system is the coordinate of Subject mapped to Query, the vertical coordinate represents the alignment coordinate of Ref1, the horizontal coordinate represents the alignment coordinate of Ref2 HSP mapped to Ref1, the matrix is constructed according to the new coordinate sorting, if two regions can be aligned, the alignment score is put into the corresponding position of the matrix; S3, for any one cell, using dynamic programming to HSP in series, the formula is as follows: ; where: i is the current HSP, j is the previous HSP, is the score of the current HSP, is the interval penalty when connecting i and j; S4. Scanning the overall matrix for the highest scoring value and using a dynamic programming algorithm to calculate the sum of base coordinates involved in the path of the highest value X ; S5、Calculate X If the ratio of Query to Subject is greater than a set threshold, Query is marked as a homologous sequence of Subject and removed. S6. For any one Contig set S, the above steps S1-S5 algorithm splits the Contig into two sets, and ; S7, for the logical loophole that appears from the self-comparison elimination method, from S 保留 any one C, query the previous dynamic programming results, find the S 保留 containing the first two Contigs C x and C y , C x and C y concatenate the sequences to obtain the sequence C connect , C connect align with C, and then perform steps S1-S6 on both and update S 保留 ; S8, traverse S 保留 Until end.

Images