Methods for detecting genomic structural variations between distant species

Abstract

The application discloses a method for detecting genomic structural variations between distant species, comprising the following steps: step S1, preprocessing the genomes of different species to be detected which meet the requirements of distant species; step S2, taking the reference species genome as a benchmark, and performing pairwise alignment by using LASTZ; step S3, integrating the obtained alignment results by using MULTIZ, and generating a multi-genome alignment task processing alignment task file in batches; step S4, generating an alignment block optimization task in batches, and performing filtering and realignment of alignment blocks on a multi-sequence alignment file; step S5, in a high-confidence MAF file, identifying insertion, deletion, inversion, intrachromosomal translocation and interchromosomal translocation operations, and generating a BED format variation result file; and step S6, a program generates a variation detection task to calculate the average consistency of the flanking sequences of each structural variation breakpoint in the variation result file. The method is suitable for efficient, accurate and low resource consumption genomic structural variation detection of distant species.

Description

Technical Field

[0001] This application relates to the field of genome detection technology, and in particular to a method for detecting genomic structural variations among distantly related species. Background Technology

[0002] Structural variations (SVs) are generally defined as genomic alterations of at least 50 base pairs (bp) in the sequence or content of a DNA sequence, primarily including insertions, deletions, duplications, inversions, and translocations. As a major source of genetic variation within and between species, structural variations play a crucial role in numerous research fields. Numerous studies have confirmed that structural variations are closely related to the occurrence and development of human diseases, influencing their pathogenesis, diagnosis, and treatment outcomes. In plant and animal breeding, structural variations serve as important molecular markers, guiding the breeding and selection of superior varieties. Simultaneously, structural variations are also key drivers of species evolution, significantly impacting adaptive evolution and phenotypic differentiation.

[0003] Therefore, accurate detection of structural variations in the genome is crucial for a deeper understanding of their role in disease development and trait expression. Developing more efficient and accurate methods for detecting structural variations is a key prerequisite for revealing their biological value.

[0004] Currently, mainstream methods for detecting structural variations are mainly divided into two categories: de novo assembly and read alignment. Among them, de novo genome assembly-based detection methods are widely considered the optimal strategy for identifying structural variations across species and have been widely applied in studies of multiple taxa. For example, the strategy of identifying lineage-specific structural variations through multi-species genome alignment has been successfully applied in studies of apes, hominids, and marine mammals, fully demonstrating that assembly-based structural variation research has clear application needs and broad application prospects.

[0005] In recent years, the widespread application of second-generation sequencing (NGS) and third-generation sequencing (TGS) technologies has made it possible to assemble individual and species genomes at low cost, providing a solid data foundation for assembly-based structural variation detection methods. Simultaneously, whole-genome alignment algorithms have also seen significant development and maturation, giving rise to a variety of efficient alignment tools. Examples include LASTZ, minimap2, and multi-genome alignment tools developed based on these, such as Minigraph, Minigraph-Cactus, Progressive Cactus, and MultiZ. These tools provide reliable algorithmic support for assembly-based structural variation detection. Driven by this technological background, assembly-based structural variation detection methods have made significant progress, with the emergence of detection tools such as SYRI, SVIM-asm, paftools.js, MUM&Co, Assemblytics, and Smartie-SV. However, these existing tools are all designed and developed for closely related species or populations within the same species. When applied to structural variation detection in distantly related species with significant genomic sequence differences (≥15%), they have obvious limitations. Specifically, the alignment software cannot generate effective collinear alignment blocks, causing the SV detection tools to fail. Furthermore, methods suitable for detecting structural variations among multiple distantly related species are extremely scarce. The find Regions Exclusively In Group program in Progressive Cactus is one of the few tools that can handle distantly related species. Although this tool is powerful, it suffers from high computational resource consumption and a high dependence on the quality of the guide tree, making it difficult to meet the practical needs for efficient and accurate detection of structural variations among distantly related species.

[0006] The information disclosed in the background section is intended only to enhance the understanding of the overall background of the invention and should not be construed as an admission or in any way implying that such information constitutes prior art known to those skilled in the art. Summary of the Invention

[0007] This application addresses the aforementioned technical problem by providing a method for detecting genomic structural variations among distantly related species. This method is applicable to the detection of genomic structural variations among distantly related species in a highly efficient, accurate, and low-resource-consumption manner.

[0008] This application provides a method for detecting genomic structural variations among distantly related species, comprising the following steps: Step S1: Call the 01_pre.sh program to preprocess the genomes of different species that meet the requirements of distantly related species. The genomes of different species include: target species genome, reference species genome, and background species genome. Generate standardized file sets of background species genome, target species genome, and reference species genome respectively. Step S2: Using the reference species genome as a benchmark, perform pairwise alignments of the background species genome with the reference species genome and the target species genome with the reference species genome using LASTZ. The alignment results are stored in the alignment task file. Step S3: Using the alignment results obtained by MULTIZ integration, call the 03_MSA.sh program to generate multiple genome alignment task files in batches, and then generate MAF format multiple sequence alignment files. Step S4: Call the 04_SV.sh program to generate alignment block optimization tasks in batches, filter and re-align the alignment blocks of the multi-sequence alignment files, improve the alignment quality, and obtain a high-confidence MAF file containing high-confidence alignment blocks; Step S5: In the high-confidence MAF file, identify insertion, deletion, inversion, intrachromosomal translocation, and interchromosomal translocation operations, and generate a variant result file in BED format; Step S6: Call the 04_SV.sh program to generate a batch of mutation detection tasks, calculate the average consistency of the flanking sequences of each structural mutation breakpoint in the mutation result file, and output a BED file containing the consistency value; Step S7: Obtain the specific structural variation dataset of the target species genome. The specific structural variation dataset includes: DEL, INS, INV, intra-TRA, and inter-TRA data of the target species genome relative to the reference species genome; and DEL, INS, INV, intra-TRA, and inter-TRA data of the background species genome relative to the reference species genome.

[0009] Preferably, the preprocessing in step S1 includes the following steps: Step S11: Record the length information of each chromosome in the genome to be tested using the faSize tool; Step S12: Generate a 2-bit file of the genome to be tested using the faToTwoBit tool; Step S13: Using the faSplit tool, the 2-bit file of the genome to be tested is split into chromosome sequences to obtain single chromosome files. At the same time, a ref_chr folder is automatically created to store the chromosome sequence file of the reference species' genome1.fa separately. Then, using the faToNib tool, each chromosome sequence file is converted into a nib format file, and a nib folder is automatically created for each species' genome to be tested to store the chromosome sequence file and its nib format file for each genome. Step S14: Check the chromosome sequence files of the reference species genome in ref_chr, manually delete sequence files that do not have chromosomes attached, and do not add suffixes to the sequence file names to obtain a standardized file set. The standardized file set includes: the sizes file, 2-bit file, nib folder, and the required reference species genome chromosome files in the ref_chr directory for each species genome.

[0010] Preferably, step S2 includes the following steps: Step S21: Call the 02_align.sh program to generate alignment tasks for each pair of genomes in the standardized file set in batches. Prepare the template file Template1_lastz_script.sh. Use the default parameters to automatically generate the shell_sh folder to store the alignment task files. Perform alignment tasks on the data in each folder separately, and compare each chromosome of the background species genome with each chromosome of the reference species genome, and each chromosome of the target species genome with each chromosome of the reference species genome. Step S22: After the pairing alignment task is completed, an axt folder is automatically generated and the pairing alignment result files are stored in the automatically generated axt folder, with one folder corresponding to each chromosome.

[0011] Preferably, step S3 includes the following steps: Step S31: Connect the alignment strands: Based on the axtChain tool, connect all the axt format alignment results corresponding to one chromosome of the reference species into a chain format alignment strand file, with the parameters set to "-minScore=5000 -linearGap=medium"; Step S32: Alignment chain networking: Based on the chainPreNet tool, redundant alignment chains are filtered out, and high-quality alignment chains are retained. Then, based on the chainNet tool, the filtered alignment chains are constructed into an alignment network file in net format. Step S33: Identify alignment networks with similar origins: Based on the net Syntenic tool, identify and retain alignment networks with similar origins, generate syntenic net files, and remove invalid alignment information without similar origins. Step S34, Format Conversion: Use the netToAxt tool to convert the obtained syntenic net file into an axt format homology alignment result, and then use the axtToMaf tool to convert the axt format file into a MAF format multiple sequence alignment file.

[0012] Preferably, step S4 includes the following steps: Step S41: Initialize filtering: Set the alignment block length threshold to 100bp, while retaining alignment blocks that exist in the genomes of all species to be tested, and generate chromosome i_filter.maf files for the alignment blocks of the genomes of the species to be tested respectively; Step S42, Re-classification and Re-comparison: For the comparison blocks retained after filtering, perform the following re-comparison operation based on the comparison block length: When the alignment block length is ≤1kbp, call the MUSCLE tool, set the default parameters, and perform multiple sequence re-alignment; When the alignment block length is greater than 1kbp, the MAFFT tool is invoked, the default parameters are set, and multiple sequence re-alignment is performed. After the classification and re-comparison are completed, a high-confidence MAF file containing high-confidence comparison blocks is generated.

[0013] Preferably, step S5 includes the following steps: Step S51: Insertion variant identification: Detect sequence continuity differences between the target species genome and the background species genome in multiple sequence alignment. When the sequence length of the target species is ≥50bp and the corresponding coordinate of the background species is a deletion or the interval distance is ≤10bp, it is determined to be an insertion variant. Record the variant location and the length of the variant gene sequence. Step S52, Deletion Variation Identification: When the sequence length of the background species genome is ≥50bp and the proportion of deletion symbols is ≤20%, and the gene in the target species genome is deleted or the interval distance is ≤10bp at the corresponding sequence position, it is determined to be a deletion variation, and the variation position and the length of the variation gene sequence are recorded. Step S53, Inversion Variation Identification: Detect the strand direction of the alignment blocks. When, in any three alignment blocks, the strand direction of the middle alignment block of the target species genome is opposite to that of the first and last alignment blocks, and the interval between alignment blocks of all species is ≤100bp, and the chromosome name and strand direction of the background species are consistent, it is determined to be an inversion variation. Record the variation location and the length of the variation gene sequence. Step S54, Translocation Variation Identification: Detect the chromosome names and positional arrangement of any three alignment blocks. When the chromosome name in the middle alignment block is different from that in the first and last alignment blocks, it is determined to be an interchromosomal translocation. When the chromosome name in the middle alignment block is the same as that in the first and last alignment blocks, but the position arrangement is abnormal, and the background species alignment block interval is ≤10bp, and the chromosome name and chain direction are consistent, it is determined to be an intrachromosomal translocation, and the location and length information of the mutation are recorded. After mutation identification is completed, a mutation result file in BED format is generated.

[0014] Preferably, step S6 includes the following steps: Step S61: If the mutation start coordinate and the termination coordinate are both located within the same alignment block, calculate the mean nucleotide sequence consistency of all species within the range upstream of the mutation start position and downstream of the termination position within the alignment block, and then calculate the overall mean consistency on both sides. Step S62: If the mutation start coordinates and termination coordinates are located in two different alignment blocks, calculate the mean nucleotide sequence consistency of all species in the alignment block where the mutation start position is located and the alignment block where the mutation termination position is located, respectively. Then calculate the overall mean consistency of the two alignment blocks and output a BED file containing the consistency values.

[0015] Another aspect of this application provides an apparatus for detecting genomic structural variations among distantly related species, as described above, comprising: The preprocessing module is used to: call the 01_pre.sh program to preprocess the genomes of different species that meet the requirements of distantly related species. The genomes of different species include: target species genome, reference species genome, and background species genome, and generate standardized file sets of background species genome, target species genome, and reference species genome respectively. The LASTZ pairwise alignment module is used to: use the reference species genome as a benchmark to perform pairwise alignments between the background species genome and the reference species genome, and between the target species genome and the reference species genome, and store the alignment results in the alignment task file. The MULTIZ multiple sequence integration module is used to: use the alignment results obtained by MULTIZ integration, call the 03_MSA.sh program to generate multiple genome alignment task files in batches, process the alignment task files, and generate MAF format multiple sequence alignment files; The high-confidence alignment module is used to: call the 04_SV.sh program to generate alignment block optimization tasks in batches, filter and re-align alignment blocks of multi-sequence alignment files, improve alignment quality, and obtain high-confidence MAF files containing high-confidence alignment blocks; The variant identification module is used to identify insertion, deletion, inversion, intrachromosomal translocation, and interchromosomal translocation operations in high-confidence MAF files and generate variant result files in BED format. The average consistency module is used to: call the 04_SV.sh program to generate batch mutation detection tasks, calculate the average consistency of the flanking sequences of each structural mutation breakpoint in the mutation result file, and output a BED file containing the consistency values.

[0016] The beneficial effects that this application can produce include: 1) High F1 score under high dissimilarity: This method achieves accurate detection of structural variations in target species among distantly related species through standardized preprocessing, paired genome alignment, multiple sequence alignment integration, alignment block optimization, structural variation identification, and consistency verification. The method has a wide effective detection range of 15%–45%, and simulation results show an F1 score of 30–80% at 35%–45% dissimilarity.

[0017] 2) Applicable to the detection of genomic structural variations among distantly related mammalian species: When used to detect human, house mouse, Norwegian rat and Eurasian red squirrel, it can detect 2063 human-specific DELs and 2488 INS. The SV signal support of INS and DEL is 100%. The DEL results have no obvious length shift, while other existing methods cannot detect specific INS. 3) Applicable to the detection of genomic structural variations among distantly related species in reptiles: When used to detect the Yangtze alligator, common gecko, European snake lizard and hot spring snake with a maximum sequence difference of 24.7%, a total of 630 hot spring snake-specific DELs and 1445 INS were detected. The SV signal support of the obtained INS and DELs was 100%, and 84.6% of the detected SV breakpoints had an offset of ≤5 bp. Other existing methods could not detect specific INS. 4) Operating costs are lower than existing methods: The peak memory usage of this method is only 3.80 GB, and the maximum virtual memory is a maximum of 15.73 GB, while existing methods require 19-23 GB of memory and have a maximum virtual memory of 60-70 GB. The CPU usage time of this method is 56-259 minutes, compared to 5284-9472 minutes by comparative tools, representing a reduction of over 95%. The CPU core count provided by this invention is stable at 19-21 cores, while the comparative tools fluctuate significantly from 33-56 cores. The actual runtime of the method provided by this invention is 1004 minutes and 2425 minutes, while existing comparative tools exceed 6792 minutes. The CPU core count used by the method provided by this invention is 57 and 28 cores, while existing tools use as many as 122 and 123 cores.

[0018] 5) Improve computational efficiency: The system corresponding method step setting function module for executing this method supports parallel or serial operation, improving computational efficiency. Attached Figure Description

[0019] Figure 1 A schematic diagram of the DistantSV-Cross framework in at least one embodiment provided in this application, wherein A represents the six main steps of the DistantSV-Cross framework; and B represents the flowchart and specific script implemented in each step of the DistantSV-Cross framework. Figure 2 The following diagrams illustrate simulated data test results in at least one embodiment provided in this application: A is a schematic diagram of insertion detection within and between alignment blocks; B is a schematic diagram of missing detection within and between alignment blocks; C is a schematic diagram of inversion detection between alignment blocks; and D is a schematic diagram of translocation detection between alignment blocks. Figure 3 The following are validation results for simulated datasets of nine different genomic sequence differences in at least one embodiment provided in this application: A represents the precision and recall of insertion processing; B represents the F1 score of insertion processing; C represents the precision and recall of deletion processing; D represents the F1 score of deletion processing; E represents the precision and recall of inversion processing; F represents the F1 score of inversion processing; G represents the precision and recall of intrachromosomal translocation processing; H represents the F1 score of intrachromosomal translocation processing; I represents the precision and recall of interchromosomal translocation processing; and J represents the F1 score of interchromosomal translocation processing. Figure 4 The figure shows the validation results of the real genome dataset of distant mammals in Experiment 2 provided in this application. A is the phylogenetic tree of the four mammal species in Experiment 2, including sequence divergence degree and divergence time; B is the structural variation detection results of each tool; C is the difference in the length of the breakpoints of the structural variations detected by each tool. Figure 5 The figure shows the validation results of the real genome dataset of distant reptiles in Experiment 3 provided in this application. A is the phylogenetic tree of the four reptile species in Experiment 3, including sequence divergence degree and divergence time; B is the structural variation detection results of each tool; C is the difference in the breakpoint length of the structural variation detected by each tool. Figure 6 The resource consumption statistics for Experiment 4 provided in this application are shown in the figure.

[0020] Figure 7 The image shows the PCR experimental verification results provided in Experiment 5 of this application. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0022] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0023] Technical means not detailed in this application and not used to solve the technical problems of this application are all set according to common general knowledge in the field, and multiple common general knowledge setting methods can be implemented.

[0024] See Figures 1-7 The method for detecting genomic structural variations among distantly related species provided in this application includes the following steps: Step S1: After calling the 01_pre.sh program to preprocess the genomes of different species that meet the requirements of distantly related species, the genomes to be detected include: the genome of the target species, the genome of the reference species, and the genome of the background species, and generate standardized file sets of the genomes of the background species, the target species, and the reference species respectively. Step S2: Using the reference species genome as a benchmark, perform pairwise alignments of the background species genome with the reference species genome and the target species genome with the reference species genome using LASTZ. The alignment results are stored in the alignment directory file. Step S3: Using the alignment catalog file obtained by MULTIZ integration, call the 03_MSA.sh program to batch construct multi-genome alignment task processing, and generate MAF format multiple sequence alignment files; Step S4: Call the 04_SV.sh program to generate alignment block optimization task files in batches, filter and re-align the alignment blocks of the multi-sequence alignment files, improve the alignment quality, and obtain MAF format multi-sequence alignment files containing high-confidence alignment blocks; Step S5: In the MAF format multiple sequence alignment file containing high-confidence alignment blocks, use modeling algorithm operations to identify insertions, deletions, inversions, intrachromosomal translocations and interchromosomal translocations, and generate a BED format variant detection result file; Step S6: Call the 04_SV.sh program to batch calculate the average consistency of the flanking sequences of each structural variant breakpoint in the variant detection result file, and output a BED file containing the consistency value; Step S7: Obtain the specific structural variation dataset of the target species genome. The specific structural variation dataset includes: DEL, INS, INV, intra-TRA, and inter-TRA data of the target species genome relative to the reference species genome; and DEL, INS, INV, intra-TRA, and inter-TRA data of the background species genome relative to the reference species genome.

[0025] Specifically, the method provided in this application includes the following steps: like Figure 1 As shown, the method provided in this application includes: sequentially performing genomic data preprocessing, pairwise genome alignment, multiple sequence alignment integration, alignment block optimization, structural variation identification, and consistency verification. Four genomes that meet the requirements of distantly related species are selected as examples (genome1.fa, genome2.fa, genome3.fa, genome4.fa), and the genomic sequence differences between the selected genomes are ≥15%. The phylogenetic relationship is "(genome1.fa genome2.fa) genome3.fa)genome4.fa)," where genome1.fa is the reference species genome, which contains 5 chromosome sequences, genome3.fa is the target species genome, and the background species genomes are genome2.fa and genome4.fa. It is worth noting that in real-world genome dataset analysis, the reference genome is selected based on the quality of the species' genome assembly and annotation. In this application's example, simulated genome data is used, and the reference genome can be either genome2.fa or genome4.fa, offering more flexibility. However, the target species' genome must be genome3.fa because it contains known structural variations compared to genome1.fa, genome2.fa, and genome4.fa. This application's example aims to detect structural variations specific to genome3.fa. The method provided in this application will then be used to detect these unique structural variations, as detailed below: Step S1: Preprocess the genomic data of multiple species to generate a standardized file set; The 01_pre.sh program is called to process the four genome files sequentially. First, the length information of each chromosome in the four genomes is recorded using the faSize tool, namely genome1.sizes, genome2.sizes, genome3.sizes, and genome4.sizes. Secondly, based on the faToTwoBit tool, four 2-bit files for the genomes are generated, namely genome1.2bit, genome2.2bit, genome3.2bit, and genome4.2bit. The efficiency of the file format is used to save storage space and facilitate fast random access to large-scale genome data. Subsequently, using the faSplit tool, the four genome files were split according to chromosome sequences to obtain individual chromosome files. In addition, a ref_chr folder was automatically created to store the chromosome sequence file of the reference species genome1.fa. Then, using the faToNib tool, each chromosome file was converted into a nib format file, and a nib folder was automatically created for each species to store the chromosome sequence file and its nib format file for each genome, namely genome1.nib, genome2.nib, genome3.nib, and genome4.nib. Finally, check the chromosome sequence files of the reference species genome in ref_chr, manually delete the sequence files that do not have chromosomes attached, and do not add a suffix to the sequence file names. In this embodiment, the contents of this folder are as follows: "chr1, chr2, chr3, chr4, chr5".

[0026] After preprocessing, a standardized file set is obtained, which includes: the sizes file, 2-bit file, nib folder, and reference species genome chromosome file required in the ref_chr directory for each species' genome. After preparation, proceed to the next step.

[0027] Step S2: Using the reference species genome as a benchmark, perform pairwise alignments using LASTZ, and store the alignment results in the alignment task file. To call the 02_align.sh program to generate alignment tasks for paired genomes in the standardized file set in batches, you first need to prepare the template file Template1_lastz_script.sh. Using the default parameters, the shell_sh folder will be automatically generated to store the alignment task files.

[0028] In this embodiment, the reference species chr1 sequence is sequentially compared with genome2.fa, genome3.fa, and genome4.fa to generate corresponding alignment tasks. The numbers in genome2.fa, genome3.fa, and genome4.fa are the genome numbers of the reference species. The remaining sequence files of the reference species are generated in the same way. Finally, 15 alignment task files (for supporting serial operation) and 1 parallel task file combined_lastz_job.sh (for supporting parallel operation) are generated in the shell_sh folder. You can flexibly choose one of the operation modes according to the actual situation.

[0029] After the alignment task is completed, an automatically generated axt folder will be created, containing the alignment result files. Each chromosome will have its own folder, namely chr1_axt, chr2_axt, chr3_axt, chr4_axt, and chr5_axt. The numbers in chr1_axt, chr2_axt, chr3_axt, chr4_axt, and chr5_axt represent the chromosome numbers in each genome. Once the alignment task files are successfully generated, you can proceed to the next step.

[0030] Step S3: Using the alignment results obtained by MULTIZ integration, call the 03_MSA.sh program to generate multi-genome alignment task files in batches, and generate multiple sequence alignment MAF files; First, a template file, Template2_multiz_script.sh, needs to be prepared, using default parameters. The alignment tasks are stored in the shell_sh folder. In this embodiment, each chromosome of the reference species corresponds to one alignment task. Ultimately, five alignment task files (for supporting serial execution) and one parallel task file, combined_multiz_job.sh (for supporting parallel execution), are generated in the shell_sh folder. You can flexibly choose one of the execution methods according to the actual situation.

[0031] The multi-genome integration process in this task includes the following steps: Step S31: Connect the alignment strands. Based on the axtChain tool, connect all the axt format alignment results corresponding to each chromosome of the reference species into a chain format alignment strand file. The parameters are set to "-minScore=5000 -linearGap=medium". Step S32: Networking of alignment chains. Based on the chainPreNet tool, redundant alignment chains are filtered out, and high-quality alignment chains are retained. Then, based on the chainNet tool, the filtered alignment chains are constructed into a net-format alignment network file with the parameter set to (stdout / dev / null). Step S33: Identify the alignment network with the same origin. Based on the netSyntenic tool, identify and retain the alignment networks with the same origin from the alignment network file generated by chainNet, generate the syntenic net file, and remove invalid alignment information with no same origin. Step S34: Format conversion. Use the netToAxt tool to convert the obtained syntenic net file into an axt format homology alignment result, and then use the axtToMaf tool to convert the axt format file into a MAF format multiple sequence alignment file.

[0032] After the multi-genome alignment task is completed, folders will be automatically generated to store the alignment results of each chromosome, namely chr1_multiz_out, chr2_multiz_out, chr3_multiz_out, chr4_multiz_out, and chr5_multiz_out. Once the multi-genome alignment files are generated successfully, you can proceed to the next step.

[0033] Step S4: Call the 04_SV.sh program to generate alignment block optimization tasks in batches. This involves filtering and re-aligning alignment blocks in MAF format multi-sequence alignment files to improve alignment quality and obtain high-confidence MAF files containing high-confidence alignment blocks. The template file Template3_SV_script.sh required for the alignment block optimization task specifically includes the following operations: Step S41: Initialize filtering: Set the comparison block length threshold to 100bp, and retain the comparison blocks that exist in all 4 species, generating chr1_filter.maf, chr2_filter.maf, chr3_filter.maf, chr4_filter.maf, and chr5_filter.maf; Step S42, Re-alignment: For the retained alignment blocks after filtering, perform the following re-alignment operations based on the block length: When the block length is ≤1kbp, call the MUSCLE tool to perform multi-sequence re-alignment; when the block length is >1kbp, call the MAFFT tool with parameters set to (--inputorder --preservecase) to perform multi-sequence re-alignment; by selecting different tools according to the block length, ensure the alignment accuracy of alignment blocks of different lengths. After completing the re-alignment, generate high-confidence MAF files containing high-confidence alignment blocks: chr1_filter_realign.maf, chr2_filter_realign.maf, chr3_filter_realign.maf, chr4_filter_realign.maf, chr5_filter_realign.maf.

[0034] Step S5: In the high-confidence MAF file, identify insertion, deletion, inversion, intrachromosomal translocation, and interchromosomal translocation operations, and generate a variant result file in BED format; The 04_SV.sh program is used to generate batch variant detection tasks. The template file Template3_SV_script.sh from step S4 is used to detect structural variants (insertions, deletions, inversions, intrachromosomal translocations, and interchromosomal translocations) specific to the target species genome3.fa. The specific identification process is as follows: Step S51, Insertion variant identification: Detect the sequence continuity differences between the target species genome3.fa and the background species genome1.fa, genome2.fa, and genome4.fa in multiple sequence alignment. When the target species sequence length is ≥50bp and the corresponding coordinate of the background species is missing or the interval distance is ≤10bp, it is determined to be an insertion variant. Record the variant location and variant sequence length. Step S52, Deletion Variation Identification: When the sequence length of the background species genome1.fa, genome2.fa, and genome4.fa is ≥50bp and the proportion of missing symbols is ≤20%, and the corresponding position of the target species genome3.fa is missing or the interval distance is ≤10bp, it is determined to be a deletion variation, and the variation position and variation sequence length are recorded. Step S53, Inversion Variation Identification: Detect the chain direction of the alignment blocks. When the chain direction of the middle alignment block of the target species genome3.fa is opposite to that of the first and last alignment blocks in any of the three alignment blocks selected, and the interval between the alignment blocks of all species is ≤100bp and the chromosome name and chain direction of the background species are consistent, it is determined to be an inversion variation. Record the variation location and variation sequence length. Step S54, Translocation Variation Identification: Detect the chromosome names and positional arrangement of any three alignment blocks. When the chromosome name in the middle alignment block is different from that in the first and last alignment blocks, it is determined to be an interchromosomal translocation. When the chromosome name in the middle alignment block is the same as that in the first and last alignment blocks, but the position arrangement is abnormal, and the background species alignment block interval is ≤10bp, and the chromosome name and chain direction are consistent, it is determined to be an intrachromosomal translocation, and information such as the location and length of the mutation is recorded. After the mutation identification module finishes running, it generates a mutation result file in BED format, which is then used as input for the next step.

[0035] Step S6: Call the 04_SV.sh program to generate a batch of mutation detection tasks, calculate the average consistency of the flanking sequences of each structural mutation breakpoint in the mutation result file, and output a BED file containing the consistency value. Based on the template file Template3_SV_script.sh mentioned above, the consistency of sequences on both sides of the mutation breakpoint is calculated for the mutation result file in BED format.

[0036] The calculation process includes the following steps: Step S61: If the mutation start coordinate and the termination coordinate are both located within the same alignment block, calculate the mean nucleotide sequence consistency of all species within the range upstream of the mutation start position and downstream of the termination position within the alignment block, and then calculate the overall mean consistency on both sides. Step S62: If the mutation start coordinates and termination coordinates are located in two different alignment blocks, calculate the mean nucleotide sequence consistency of all species in the alignment block where the mutation start position is located and the alignment block where the mutation termination position is located, and then calculate the overall mean consistency of the two alignment blocks.

[0037] In the above six steps of this invention, steps S1, S2, and S3 are executed independently, and after being executed sequentially, they proceed to steps S4-6, which run in the same task. In this embodiment, five tasks are generated in the shell_sh file, each representing one of the five chromosomes of the reference species genome. Each task integrates steps S4-6, and includes five serially executed tasks and one parallelly executed task. Finally, the calculated mean nucleotide sequence identity is added to a BED format variant file, which contains complete variant information.

[0038] Step S7: Sample detection results and obtain target species-specific structural variation dataset After performing the above six steps sequentially, a specific structural variation dataset of the target species' genome3.fa will be obtained. This dataset is located in the multiz_out subfiles of five folders, representing the five chromosomes of the reference species' genome: "chr1_multiz_out / multiz_out, chr2_ multiz_out / multiz_out, chr3_multiz_out / multiz_out, chr4_ multiz_out / multiz_out, chr5_multiz_out / multiz_out". Each folder contains five variation files, ending with the suffixes "_ins.bed, _del.bed, _inv.bed, _intraTRA.bed, _interTRA.bed". These are the final BED format target files, which record detailed information such as the location, strand direction, variation length, and upstream and downstream sequence consistency of the variation breakpoints in the reference and target species genomes. Subsequent analysis can then be conducted based on this file.

[0039] Simulation Examples The following simulation tests evaluate and compare the method provided in this application (DistantSV-Cross) with seven existing methods for detecting genomic structural variations from five aspects. The existing methods are findRegionsExclusivelyInGroup, SYRI, MUM&Co, Smartie-SV, Assemblytics, SVIM-asm, and paftools.js. The simulation results show that, compared with existing detection technologies, this invention has significant advantages in detecting structural variations among distantly related species.

[0040] Experiment 1: Validation results using simulated datasets of 9 different genomic sequence dissimilarity levels. Nine simulated datasets of genomic sequence dissimilarity levels (15%, 20%, 25%, 30%, 35%, 40%, 45%, 50%, 55%) were used. This represents a gradient of distant genomic sequence dissimilarity levels, which can be used to evaluate the ability of this invention and existing technologies to detect structural variations in distant dissimilarity gradients.

[0041] The obtained test results are as follows Figure 3As shown, DistantSV-Cross exhibits the best overall performance and the widest applicability across sequence divergences. Specifically, in terms of precision and recall, within the 15%-30% sequence divergence range, DistantSV-Cross demonstrates significantly better precision (73.65%~88.89%) and recall (95.24%~97.89%) for deletions (87.20%~97.60%) and insertions (80.00%~92.80%) than other software. Within the 15%-45% range, the precision for inversions, intrachromosomal translocations, and interchromosomal translocations reaches 99.69%~100%, failing only when the sequence divergence is ≥50%. In comparison, SYRI, MUM&Co, Smartie-SV, and Assemblytics tools failed to detect any structural variations when the divergence was ≥15%; SVIM-asm and paftools.js could only operate effectively at 15% and 20% with low precision and recall, respectively, and became completely ineffective at ≥25%; findRegionsExclusivelyInGroup had a high missing recall (around 85%), but poor precision for missing (9%~21%) and insertion (3%~5%) data.

[0042] Regarding F1 scores, DistantSV-Cross's overall F1 score gradually decreased from 91.98% at 15% to 32.06% at 45%, making it the only tool that consistently performed effectively between 15% and 45%. Other tools either had no F1 score after 15% or only a few limited F1 scores at low divergence levels; findRegionsExclusivelyInGroup exhibited low precision across all gradients, resulting in poor F1 performance.

[0043] In summary, the test results show that the effective detection range of the DistantSV-Cross method provided by this invention is 15% to 45%, and it is basically undetectable for 50% and above. The DistantSV-Cross method provided by this invention is significantly superior to other comparative tools.

[0044] Experiment 2: Validation results using a real mammalian genome dataset. Four mammalian species—human, mouse, rat, and squirrel—were selected, with a maximum sequence difference of 19.1% (corresponding to approximately 87 million years of divergence), to evaluate the ability to detect structural variations in real, distantly related genomes. The results are as follows: Figure 4 As shown: Only the method provided by this invention and the findRegionsExclusivelyInGroup tool successfully detected the changes; the other six tools (SYRI, SVIM-asm, paftools.js, MUM&Co, Assemblytics, and Smartie-SV) failed to detect any valid structural variations.

[0045] The method provided by this invention detected 2063 human-specific deletions and 2488 insertions, while findRegionsExclusivelyInGroup detected only 4661 deletions and no insertions. Further statistical results show that the SV signal support for both insertions and deletions in this invention is 100%, while findRegionsExclusivelyInGroup has no insertions available for SV signal support checking.

[0046] Furthermore, this invention significantly outperforms other methods in terms of length offset, while the missing data detected by findRegionsExclusivelyInGroup exhibits a large number of length offsets, indicating that this method has a significant error in determining the missing length.

[0047] In summary, in real mammalian datasets, this invention demonstrates significant advantages in both SV signal support and length offset.

[0048] Experiment 3: Validation results using a real reptile genome dataset. Four reptile species—the Chinese alligator, the common gecko, the European snake lizard, and the hot spring snake—were selected, with a maximum sequence difference of 24.7% (corresponding to approximately 280 million years of divergence)—to evaluate the ability to detect real distant genome structural variations at higher sequence difference levels.

[0049] The results are as follows Figure 5 As shown: Only the method provided by this invention and the findRegionsExclusivelyInGroup tool successfully detected the variants; the other six tools (SYRI, SVIM-asm, paftools.js, MUM&Co, Assemblytics, and Smartie-SV) failed to detect any valid structural variations. This invention detected 630 deletions and 1445 insertions specific to *Sinocyclocheilus spp.*, while findRegionsExclusivelyInGroup detected only 81 deletions and no insertions. Further statistical results show that the SV signal support for both insertions and deletions in this invention is 100%, while findRegionsExclusivelyInGroup has no insertions available for checking SV signal support. Furthermore, this invention significantly outperforms in length offset performance, with 84.6% of detected SV breakpoint offsets ≤5bp; while the deletions detected by findRegionsExclusivelyInGroup exhibit a large number of length offsets, with only 18.5% reaching this threshold, indicating a significant error in the method's deletion length determination.

[0050] In summary, in real crawling datasets, this invention also demonstrates significant advantages in SV signal support and length offset, and its detection capability is more stable under high dissimilarity.

[0051] Experiment 4: Comparison and verification results of computational resource consumption for various technologies. For example... Figure 6 As shown, on simulated datasets (sequence difference 15%~45%) and real datasets (mammals and crawlers), the five key computational metrics (maximum physical memory, maximum virtual memory, number of CPU cores, CPU time, and actual runtime) of the method provided in this invention and findRegionsExclusivelyInGroup are evaluated.

[0052] The results show that the resource consumption of the present invention is significantly lower than that of the comparison tool on all datasets. In the simulated data, the maximum physical memory of the present invention is consistently below 4 GB (peak 3.80 GB), while the comparison tool requires 19-23 GB; the maximum virtual memory of the present invention is at most 15.73 GB, less than a quarter of that of the comparison tool (60-70 GB); the CPU usage time of the method provided by the present invention is 56-259 minutes, while that of the comparison tool is 5284-9472 minutes, a reduction of more than 95%; the actual running time of the method provided by the present invention is 25-55 minutes, while that of the comparison tool is 152-180 minutes; the number of CPU cores provided by the present invention is stable at 19-21 cores, while that of the comparison tool fluctuates significantly from 33-56 cores.

[0053] In real mammalian and crawler datasets, this invention also demonstrates a significant advantage: the maximum physical memory required is 79.0 GB and 38.9 GB, respectively, while existing tools exceed 177 GB. The actual running time of the method provided by this invention is 1004 minutes and 2425 minutes, while existing comparison tools exceed 6792 minutes; the number of CPU cores used in the method provided by this invention is 57 cores and 28 cores, while existing tools have as many as 122 cores and 123 cores.

[0054] In summary, this invention demonstrates superior computational efficiency across all five metrics, consumes fewer resources, and runs faster, making it particularly suitable for large-scale analyses involving multiple distantly related species.

[0055] Experiment 5: Validation results of the PCR experiment.

[0056] Based on the detection results of Experiment 2, namely the validation results of the real mammalian genome dataset, a total of 2063 human-specific deletions and 2488 human-specific insertions were obtained. From these, 20 deletions (8 deletions in exons, 6 deletions in introns, 4 deletions in the upstream 2kb region of the gene, and 2 deletions in the distant intergenic region) and 20 insertions (8 insertions in exons, 6 insertions in introns, 4 insertions in the upstream 2kb region of the gene, and 2 insertions in the distant intergenic region) were selected for PCR experiments and Sanger sequencing. The specific operational steps are as follows: Genomic DNA was extracted from human HaCaT cells, mouse, rat, and red squirrel muscle tissue using the TIANamp kit as the DNA template for PCR experiments. Four sets of primers were designed for each variant across four species using Primer Premier 5. After PCR amplification confirming that the band length matched the variant, the product was purified and Sanger paired-end sequencing was performed on an ABI 3730XL platform. The sequencing results were aligned using MEGA11, and the alignment results were used as the true values ​​to evaluate the accuracy of structural variant detection in this invention. The obtained PCR verification results are as follows: Figure 7 As shown, after Sanger sequencing, all 40 structural variants were verified as human-specific structural variant sequences, indicating that the specific gene source obtained by the method provided in this application is accurate, and further demonstrating the reliability of the variant detection results of this invention.

[0057] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for detecting genomic structural variations among distantly related species, characterized in that, Includes the following steps: Step S1: Call the 01_pre.sh program to preprocess the genomes of different species that meet the requirements of distantly related species. The genomes of different species include: target species genome, reference species genome, and background species genome. Generate standardized file sets of background species genome, target species genome, and reference species genome respectively. Step S2: Using the reference species genome as a benchmark, perform pairwise alignments of the background species genome with the reference species genome and the target species genome with the reference species genome using LASTZ. The alignment results are stored in the alignment task file. Step S3: Using the alignment results obtained by MULTIZ integration, call the 03_MSA.sh program to generate multiple genome alignment task files in batches, and then generate MAF format multiple sequence alignment files. Step S4: Call the 04_SV.sh program to generate alignment block optimization tasks in batches, filter and re-align the alignment blocks of the multi-sequence alignment files, improve the alignment quality, and obtain a high-confidence MAF file containing high-confidence alignment blocks; Step S5: In the high-confidence MAF file, identify insertion, deletion, inversion, intrachromosomal translocation, and interchromosomal translocation operations, and generate a variant result file in BED format; Step S6: Call the 04_SV.sh program to generate a batch of mutation detection tasks, calculate the average consistency of the flanking sequences of each structural mutation breakpoint in the mutation result file, and output a BED file containing the consistency value; Step S7: Obtain the specific structural variation dataset of the target species genome. The specific structural variation dataset includes: DEL, INS, INV, intra-TRA, and inter-TRA data of the target species genome relative to the reference species genome; and DEL, INS, INV, intra-TRA, and inter-TRA data of the background species genome relative to the reference species genome.

2. The method for detecting genomic structural variations among distantly related species according to claim 1, characterized in that, The preprocessing in step S1 includes the following steps: Step S11: Record the length information of each chromosome in the genome to be tested using the faSize tool; Step S12: Generate a 2-bit file of the genome to be tested using the faToTwoBit tool; Step S13: Using the faSplit tool, the 2-bit file of the genome to be tested is split into chromosome sequences to obtain single chromosome files. At the same time, a ref_chr folder is automatically created to store the chromosome sequence file of the reference species' genome1.fa separately. Then, using the faToNib tool, each chromosome sequence file is converted into a nib format file, and a nib folder is automatically created for each species' genome to be tested to store the chromosome sequence file and its nib format file for each genome. Step S14: Check the chromosome sequence files of the reference species genome in ref_chr, manually delete sequence files that do not have chromosomes attached, and do not add suffixes to the sequence file names to obtain a standardized file set. The standardized file set includes: the sizes file, 2-bit file, nib folder, and the required reference species genome chromosome files in the ref_chr directory for each species genome.

3. The method for detecting genomic structural variations among distantly related species according to claim 1, characterized in that, Step S2 includes the following steps: Step S21: Call the 02_align.sh program to generate alignment tasks for each pair of genomes in the standardized file set in batches. Prepare the template file Template1_lastz_script.sh. Use the default parameters to automatically generate the shell_sh folder to store the alignment task files. Perform alignment tasks on the data in each folder separately, and compare each chromosome of the background species genome with each chromosome of the reference species genome, and each chromosome of the target species genome with each chromosome of the reference species genome. Step S22: After the pairing alignment task is completed, an axt folder is automatically generated and the pairing alignment result files are stored in the automatically generated axt folder, with one folder corresponding to each chromosome.

4. The method for detecting genomic structural variations among distantly related species according to claim 1, characterized in that, Step S3 includes the following steps: Step S31: Connect the alignment strands: Based on the axtChain tool, connect all the axt format alignment results corresponding to one chromosome of the reference species into a chain format alignment strand file, with the parameters set to "-minScore=5000 -linearGap=medium"; Step S32: Alignment chain networking: Based on the chainPreNet tool, redundant alignment chains are filtered out, and high-quality alignment chains are retained. Then, based on the chainNet tool, the filtered alignment chains are constructed into an alignment network file in net format. Step S33: Identify alignment networks with similar origins: Based on the net Syntenic tool, identify and retain alignment networks with similar origins, generate syntenic net files, and remove invalid alignment information without similar origins. Step S34, Format Conversion: Use the netToAxt tool to convert the obtained syntenic net file into an axt format homology alignment result, and then use the axtToMaf tool to convert the axt format file into a MAF format multiple sequence alignment file.

5. The method for detecting genomic structural variations among distantly related species according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Initialize filtering: Set the alignment block length threshold to 100bp, while retaining alignment blocks that exist in the genomes of all species to be tested, and generate chromosome i_filter.maf files for the alignment blocks of the genomes of the species to be tested respectively; Step S42, Re-classification and Re-comparison: For the comparison blocks retained after filtering, perform the following re-comparison operation based on the comparison block length: When the alignment block length is ≤1kbp, call the MUSCLE tool, set the default parameters, and perform multiple sequence re-alignment; When the alignment block length is greater than 1kbp, the MAFFT tool is invoked, the default parameters are set, and multiple sequence re-alignment is performed. After the classification and re-comparison are completed, a high-confidence MAF file containing high-confidence comparison blocks is generated.

6. The method for detecting genomic structural variations among distantly related species according to claim 1, characterized in that, Step S5 includes the following steps: Step S51: Insertion variant identification: Detect sequence continuity differences between the target species genome and the background species genome in multiple sequence alignment. When the sequence length of the target species is ≥50bp and the corresponding coordinate of the background species is a deletion or the interval distance is ≤10bp, it is determined to be an insertion variant. Record the variant location and the length of the variant gene sequence. Step S52, Deletion Variation Identification: When the sequence length of the background species genome is ≥50bp and the proportion of deletion symbols is ≤20%, and the gene in the target species genome is deleted or the interval distance is ≤10bp at the corresponding sequence position, it is determined to be a deletion variation, and the variation position and the length of the variation gene sequence are recorded. Step S53, Inversion Variation Identification: Detect the strand direction of the alignment blocks. When, in any three alignment blocks, the strand direction of the middle alignment block of the target species genome is opposite to that of the first and last alignment blocks, and the interval between alignment blocks of all species is ≤100bp, and the chromosome name and strand direction of the background species are consistent, it is determined to be an inversion variation. Record the variation location and the length of the variation gene sequence. Step S54, Translocation Variation Identification: Detect the chromosome names and positional arrangement of any three alignment blocks. When the chromosome name in the middle alignment block is different from that in the first and last alignment blocks, it is determined to be an interchromosomal translocation. When the chromosome name in the middle alignment block is the same as that in the first and last alignment blocks, but the position arrangement is abnormal, and the background species alignment block interval is ≤10bp, and the chromosome name and chain direction are consistent, it is determined to be an intrachromosomal translocation, and the location and length information of the mutation are recorded. After mutation identification is completed, a mutation result file in BED format is generated.

7. The method for detecting genomic structural variations among distantly related species according to claim 1, characterized in that, Step S6 includes the following steps: Step S61: If the mutation start coordinate and the termination coordinate are both located within the same alignment block, calculate the mean nucleotide sequence consistency of all species within the range upstream of the mutation start position and downstream of the termination position within the alignment block, and then calculate the overall mean consistency on both sides. Step S62: If the mutation start coordinates and termination coordinates are located in two different alignment blocks, calculate the mean nucleotide sequence consistency of all species in the alignment block where the mutation start position is located and the alignment block where the mutation termination position is located, respectively. Then calculate the overall mean consistency of the two alignment blocks and output a BED file containing the consistency values.

Images