With the whole genome data of Salix sinopurpurea (Spu) and Brassica rapa (Bra) as the main research materials, comparative genomics was used to compare the collinearity between Salix sinopurpurea and the reference genome Vitis vinifera (Vvi), Brassica rapa and the reference genome Arabidopsis thaliana (Ath).

Genomes and their gene annotations of Salix sinopurpurea and Vitis vinifera were downloaded from Joint Genome Institute. Download the required documents for Brassica rapa and Arabidopsis thaliana at http://brassicadb.org/and http://www.arabidopsis.org/shangxia, respectively.

Due to the huge amount of original genome data, in order to extract target data from the genome sequence and annotation files, the downloaded genome data is processed with a custom python script to obtain the blast results, which is convenient for subsequent research and analysis. Screen the original data of species genomes, information was extracted from the genome annotation files, which include chromosome number, gene start and end positions, gene transcription direction, and gene ID information, and then rename the gene ID and the number of the genes was given in order of their appearance on chromosomes. Map the gene ID in the CDS sequence and protein sequence file to the new ID of the corresponding gene in the genome annotation file. Label the processed genomic data with a unified naming method.

Blastp was used to explore to align genomic sequences of different species. Screen out gene pairs with the expected value (E-value) not greater than 10–5 and score evaluation (Score) higher than 100, so that the subsequent genome collinearity analysis results are more reliable.

The WGDI (Sun et al., 2021) use MAFFT (Wong, Suchard, and Huelsenbeck 2008) or MUSCLE (Edgar 2004) to perform multiple sequence alignment, and calculates the synonymous substitution rate using the yn00 (Yang et al., 2000) or ng86 (Nei and Gojobori 1986) program of the PAML package. Finally, the visualization is realized by extracting block, and then output blockinfo file.

In this paper, the input for DBSCAN requires the blockinfo file generated by WGDI and the chromosome length information (len file) of the two species. By setting the epsilon (eps) and minimum points (MinPts), cluster analysis is performed on the collinearity fragments in the K _S dotplot. The collinear region was then obtained from the clustering results combined with certain rules for merging. And then realizes the automatic identification and labeling of the collinear region. The comparison result of a chromosome of the target species and a chromosome of the reference species is shown as a cell on the K _S dotplot, that is, a comparison unit.

The K _S dotplot between Salix sinopurpurea and Vitis vinifera genome homologous genes drawn by wgdi (Figure 1A). The horizontal axis represents the chromosome of the target species (Salix sinopurpurea) and the vertical represents the chromosome of the reference species (Vitis vinifera). On the K _S dotplot, the chromosome number of Salix sinopurpurea is shown from left to right, and the chromosome number of Vitis vinifera is shown from top to bottom. The K _S value ranges from 0.00 to 2.00. As shown in the figure, different colored points correspond to different K _S values. It can be observed in Figure 1A, in addition to the clear and complete homologous fragments of grape chromosome 4 with Salix sinopurpurea chromosomes 6 and 18, it also has fuzzy and unclear homologous fragments with Salix sinopurpurea chromosomes 1, 2 and 4. The reason why these fragments are unclear and incomplete is that they are doubled by the whole genome triplication events shared by older dicots. The collinearity of the homologous fragments produced by the whole genome triplication events shared by ancient dicotyledons is far inferior to that of the whole genome duplication events shared by the Salicaceae. The specific manifestation is that the K _S value is significantly large, belongs to the blue-purple system, scarce and fragmented seriously. The results showed the structural similarities and differences between genomes. The generated data and pictures provide references for follow-up research.

Using the DBSCAN algorithm, by setting the eps to 50 and the MinPts to 3, cluster the K _S dotplot between Salix sinopurpurea and Vitis vinifera genome homologous genes. The algorithm outputs the clustering result figure (Figure 1B), in which each category is represented by a rectangular box. The DBSCAN can cluster out complete collinearity fragments, in grape chromosome 14 and Salix sinopurpurea chromosome 5, as well as in grape chromosome 14 and Salix sinopurpurea chromosome 15. It also can identify fragmented collinearity fragments in grape chromosome 4 and Salix sinopurpurea chromosome 5, grape chromosome 13 and Salix sinopurpurea chromosome 16. These will accurately reflect the relationship between the collinearity fragments and improve the subsequent combination effect.

The model sorts the category in the same comparison unit from top to bottom to generate ID, and calculates the number of homologous gene points in the box (num), the length of the box (y1-y2), and the width of the box (x1-x2). The model then generates cluster.csv files that contain the target species chromosome number (chr1), the reference species chromosome number (chr2), ID, and the horizontal and vertical coordinates of the upper left corner point are l_x and l_y, respectively, num, the horizontal and vertical coordinates of the lower right point are r_x and r_y, respectively, y1-y2, and x1-x2. Part of the data in the cluster.csv file of grape chromosome 13 is shown in Table 1. Among them, chromosome 13 and chromosome 1 form a class. The coordinates of the upper left corner of this class are 145, 354, and the coordinates of the lower right corner are 257, 228. The number of homologous genes contained is 32. The length of the cluster box is 126 coordinate lengths, and the width is 112 coordinate lengths.

In order to perform a combined analysis on the identified collinearity fragments, the model read the cluster.csv file generated by clustering. The model uses y_gap and x_gap, which represents the gap in the longitudinal and horizontal directions of adjacent collinear segments, as the basis for judging overlap. The location information of the gene is combined to set the parameters gap and Slen. In the comparison unit, the parameter gap represents the mean value of y_gap. Through a series of experiments and continuous optimization of parameter selection, it is finally determined that 1/6 of the corresponding chromosome length of the target species is the value of parameter Slen. For all collinearity fragments whose num is greater than the specified value, one condition is that the collinearity fragments do not overlap, another is overlap. In the first condition, collinear fragments will be merged if 0 ≤y_gap ≤ gap and 0 ≤x_gap ≤ Slen. And in the second condition, there is overlap in the y-axis direction, they will be merged when it meets 0 ≤ x_gap < Slen; if there is overlap in the x-axis direction, when it meets 0 ≤ y_gap ≤ gap, merge them. Finally, the model output the clustering result graph (Figure 1C), in which the merged result is marked with a rectangular box, and the combine.csv file is generated.

The content of the combine.csv file is the same as the cluster.csv file. Table 2 shows some data of the combine.csv file. Among them, chromosome 13 and chromosome 8 form two classes. The coordinates of the upper left corner of the first class are 592, 1,264, and the lower right corner are 671, 1,155. The number of homologous genes contained in this class is 56. The length of the labeled box ∆y is 112 coordinate length, and the width ∆x is 79 coordinate length; the upper left corner coordinate of the second class are 78, 1,089, and the lower right corner are 605, 446. The number of homologous genes contained in this class is 239, ∆y is 643 coordinate length, ∆x is 537.

In order to determine the polyploid types and related chromosomes of the species, we develop the polyploidy recognition algorithm. The algorithm read the generated combine.csv file, look for the labeled box with the largest ∆y, mark it, and then look up and down to find the labeled box with the length less than it. The result.csv file will be built by adopting exhaustive above process. Among them, comro represents the combination round, sumy represents the sum of ∆y in same combination round. In order to determine the specific polyploidy of the species, the gene collinearity evaluation index line chart is drawn. The horizontal is the combined round, and the vertical is the corresponding gene collinearity evaluation index. The significant inflection point in the line chart represents the corresponding polyploid type. The gene collinearity evaluation index (MI) was calculated by dividing the cumulative collinearity fragments length to the corresponding chromosome length of the reference species in the len file to describe the performance of the polyploidy in the corresponding combination round, and the larger its value, the better the performance. M I = ∑ m = 1 n Δ y m l e n i Where Δ y m and n are the length and number of collinearity fragments in the same combination round, respectively;   l e n i is the corresponding chromosome length of the reference species, i is the corresponding chromosome number.

After determining the combination round, output the combined result graph. To describe the analysis of input, output and fetching the final results of the analysis, we made pseudo code. The pseudo code of the chromosome collinearity fragment labeling and polyploidy recognition algorithm is shown as below for a better understanding of the context and better assess the relevance of this paper.

Algorithm 1

Chromosome collinearity fragment labeling and polyploidy recognition algorithm.

Using PolyReco to objectively determine the polyploid types of Salix sinopurpurea. The model use Vitis vinifera as the reference genome to identify the target species Salix sinopurpurea polyploid type, and read the data of Vitis vinifera chromosome 4, 13 and 14. The DBSCAN algorithm obtains the collinearity fragments, and get the combine.csv file. Then using the polyploidy recognition algorithm to combine the labeled boxes exhaustively, and get the gene collinearity evaluation index table of each chromosome in different combination rounds (Supplementary Table S1). In the colinearity evaluation index line chart (Figure 2A), we can find that chromosomes 4, 13, and 14 of Vitis vinifera have obvious inflection points when the combined round is 2. Therefore, it is determined that the Salix sinopurpurea has a whole genome duplication event recently. This conclusion can be found in Wang et al., 2011.

After determining the specific polyploidy, we can output the information of the labeled box participating in the polyploidy, and obtain the result.csv file of Vitis vinifera chromosomes 4, 13, and 14, in which the data of chromosome 13 is shown in Table 3.

According to the result.csv file, output the combined figure of the Salix sinopurpurea polyploidy (Figure 2B). When the combination round is 2, get the two groups with the highest scores, among which the chromosomes 5 and 6 of the Salix sinopurpurea can be combined into a relatively complete chromosome 4 of Vitis vinifera, the corresponding MI is 67.87%; the chromosomes 2 and 18 of the Salix sinopurpurea can be combined into a relatively complete chromosome 4 of Vitis vinifera, and the corresponding MI is 42.19%. The chromosomes 8 and 16 of the Salix sinopurpurea can be combined into a relatively complete Vitis vinifera chromosome 13 with MI of 73.93%; the chromosomes 9 and 10 of the Salix sinopurpurea can be combined into a relatively complete Vitis vinifera chromosome 4, MI is 70.96%. The chromosomes 5, 15, and 17 of Salix sinopurpurea can be combined into a relatively complete Vitis vinifera chromosome 14 with MI of 88.74%; the chromosomes 10, 13, 16, and 17 of Salix sinopurpurea can be combined into a relatively complete Vitis vinifera chromosome 14, MI is 79.88%.

In order to further verify the universality of the method, according to the procedure in 3.1, the model uses Arabidopsis thaliana as the reference genome to identify the polyploidy type of the target species Brassica rapa. Through reading the data of chromosomes 1, 3 and 5 of Arabidopsis thaliana, we finally get the gene collinearity evaluation index table of each chromosome in different combination rounds (Supplementary Table S2). In the collinearity evaluation index line chart (Figure 3A), it can be found that chromosomes 1, 3, and 5 of Arabidopsis thaliana had an obvious turning point when the combination round was 3. Therefore, it is determined that the Brassica rapa. had a whole genome triplication event recently. This conclusion can be found in Wang (2011a).

After determining the specific polyploidy, the result.csv file of Arabidopsis thaliana chromosomes 1, 3, and 5 is obtained, in which the data of chromosome 3 is shown in Table 4.

According to the result.csv file, output the combined figure of Brassica rapa polyploidy (Figure 3B). When the combination round is 3, get the three groups with the highest scores, among which the chromosomes 7, 8, and 9 of the Brassica rapa can be combined into a relatively complete chromosome 1 of Arabidopsis thaliana, the corresponding MI is 89.88%; the chromosomes 2, 6, and 9 of the Brassica rapa can be combined into a relatively complete chromosome 1 of Arabidopsis thaliana, and the corresponding MI is 58.31%; the chromosomes 6, 7, and 10 of the Brassica rapa can be combined into a relatively complete chromosome 1 of Arabidopsis thaliana, and the corresponding MI is 49.38%. The chromosomes 3, 6, and 9 of the Brassica rapa can be combined into a relatively complete Arabidopsis thaliana chromosome 3 with MI of 88.16%; the chromosomes 1 and 4 of the Brassica rapa can be combined into a relatively complete Arabidopsis thaliana chromosome 3 with MI of 69.87%; the chromosomes 5 and 7 of the Brassica rapa can be combined into a relatively complete Arabidopsis thaliana chromosome 3 with MI of 64.98%. The chromosomes 2, 6, and 9 of the Brassica rapa can be combined into a relatively complete Arabidopsis thaliana chromosome 5, MI is 83.54%; the chromosomes 6 and 10 of the Brassica rapa can be combined into a relatively complete Arabidopsis thaliana chromosome 5, MI is 59.84%; the chromosomes 3 of the Brassica rapa can be combined into a relatively complete Arabidopsis thaliana chromosome 5, MI is 51.17%.

In Figure 3B, there are two seemingly identical collinearity fragments among the four fragments formed by Arabidopsis thaliana chromosome 3 and Brassica rapa chromosome 4, with the naked eye. But only one is labeled and used, because it has a small number of homologous genes. So this method can break through the limitations of the human eye, and find chromosome fragments with strong collinearity, as well as provide a basis for objective judgment of polyploidy.