A method and apparatus for correcting long read assembly results using short read sequences
By dividing the long read and long assembly results into small modules and comparing and correcting short read and long sequences on small memory nodes, the problem of high equipment requirements and time-consuming in the existing technology is solved, and efficient error correction effect is achieved.
Patent Information
- Application Number
- CN202110452532.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2021-04-26
- Publication Date
- 2025-06-17
- Estimated Expiration
- 2041-04-26
AI Technical Summary
The existing methods of using short-read long sequence error correction, long-read long assembly results have high requirements for equipment, long-term and low efficiency, especially when there are no large memory nodes or the need to rent expensive large-memory cloud servers.
By dividing the long read and long assembly results into smaller modules, and comparing and correcting short read and long sequences on the divided modules, small memory nodes are used for error correction, reducing memory peaks and time-consuming.
It realizes the error correction of the long read and long assembly results of short read and long sequences on small memory nodes, reducing memory peaks and time-consuming, reducing error correction costs, and improving efficiency.
Smart Images

Figure CN113012758B_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of gene sequencing technology, and in particular to a method and device for correcting the results of long-read assembly using short-read sequences. Background Art
[0002] Due to the characteristics of ultra-long reads, the average read length of long-read PacBio and Nanopore sequencing technologies is generally 10k - 30k, which can perform high-level assembly on complex genomic regions such as highly repetitive sequences, transposon regions, and highly variable regions. Generally, the Contig (contig) N50 of the assembly results of animals and plants can reach more than 1MB. Currently, the main long-read assembly software includes PBCR, Falcon, MECAT, CANU, HGAP, etc., and these software all include self-correction and assembly functions. Since the average error rate of long-read sequences is generally 10 - 20%, these software first perform sequence self-correction, and then use the corrected sequences for assembly to obtain the assembly results. Since there may be certain single-base errors or structural variations in the assembly results, the original long-read sequences are subsequently used to correct the assembly results, and high-accuracy short-read sequences with an average error rate of less than 1% are used through Pilon software for correction to obtain the final long-read assembly results. The main process of long-read assembly is as Figure 1 shown.
[0003] When short-read sequences correct the assembly results using Pilon, Pilon software corrects the entire assembly result, resulting in high memory peak and long time consumption. Because the cost of large memory nodes above 100G is relatively high, and at the same time, many universities and companies do not often use them, do not configure large memory nodes, or there are tasks running on the large memory nodes and part of the memory has been occupied, and they need to wait for the tasks to run to completion before they can be executed, which will all affect the execution of the project. Currently, Alibaba Cloud and Huawei Cloud can provide rental services for large memory (100Gb) computers, but their prices are several times to dozens of times that of small memory (below 64Gb) nodes, which is very expensive.
[0004] The program execution script for correcting the long-read assembly result racon.fasta using short-read sequences NGS_1.fq and NGS_2.fq is as follows:
[0005] (1) bwa index racon.fasta, that is, Bwa index*.fasta builds a library for the assembly result
[0006] (2) bwa mem -t 64 -M racon.fasta NGS_1.fq NGS_2.fq > ONTmin_ref.sam
[0007] (3) samtools view - bs ONTmin_ref.sam > ONTmin_ref.bam, that is, convert the *.sam format to the *.bam format
[0008] (4) samtools sort ONTmin_ref.bam ONTmin_ref, that is, the sorting command for the *.bam file. Here, ONTmin_ref.bam is the input file, and ONTmin_ref is the prefix of the output file. The output file will finally have the suffix ".bam". That is, if ONTmin_ref is changed to A, the output result will be A.bam
[0009] (5) samtools index ONTmin_ref.bam, that is, the library building command samtools index *.bam
[0010] (6) java - Xmx300G - jar pilon - 1.22.jar -- genome racon.fasta -- ouput Pilon – outdir. / -- tracks – diploid – fix snps,indels -- threads 64 -- flank 0 -- frags ONTmin_ref.bam
[0011] Steps (1), (2), (3), (4), (5), and (6) in the program execution script respectively correspond to Figure 2 library building 101 for the long - length assembly result in Figure 2 short read sequence alignment to the long - length assembly result 102 in
[0012] conversion of the comparison result format 103 in sorting of the alignment result 104 in
[0013] library building of the alignment result 105 in
[0014] error correction of the long - length assembly result 106 in
[0015] When the assembly result is greater than 2G for a large genome, step (6) in the program execution script, that is, Figure 2 the error correction step 106 of the long - length assembly result in
[0012] requires setting 300G of memory and using 64 threads for the error correction of the entire assembly result. Its memory peak value will be very high and the time consumption will also be relatively long. Therefore, the existing method of using short read sequences to correct the long read assembly result not only has high requirements for equipment, but also is time - consuming and inefficient. Summary of the Invention
[0013] The purpose of this application is to provide an improved method and device for using short read sequences to correct the long read assembly result.
[0014] To achieve the above purpose, this application adopts the following technical solutions:
[0015] One aspect of the present application discloses a method for correcting long-read assembly results using short-read sequences, comprising the following steps:
[0016] A cutting step, including splitting the long-read assembly result at a fixed threshold to obtain the split long-read assembly result; wherein, the fixed threshold is generally set to 100 Mb. By using this fixed threshold for splitting, according to the method of the present application, only 50 G of machine memory is required to complete error correction. For example, in Experiment 2 of the embodiments of the present application, the peak memory is only 26 G, while without using the method of the present application, the peak memory is more than 200 G. In principle, the larger the fixed threshold, the larger the split, and the higher the peak memory, and vice versa. It depends on the size of the machine memory for running the method task of the present application;
[0017] A step of building a library for the long-read assembly result, including building a library for the split long-read assembly result to obtain a library of the split long-read assembly result;
[0018] A comparison step, including aligning the short-read sequences to the library of the split long-read assembly result to obtain a comparison result;
[0019] A step of matching the comparison result, including matching the comparison result obtained in the comparison step to the library of the split long-read assembly result to obtain a comparison result containing matching information;
[0020] A step of converting and sorting the comparison result, including converting the comparison result containing matching information into the bam format and sorting it; In one implementation manner of the present application, it is specifically for the racon software. Among them, the steps use bwa alignment and samtools view -bs*sam to convert to the *bam format, and the sorting uses the samtools sort command. It can be understood that if it is for other software, it can also be converted into the corresponding format, not necessarily the bam format, and no specific limitation is made here;
[0021] A step of building a library for the comparison result, including building a library for the sorted comparison result converted into the bam format to obtain a library of the comparison result;
[0022] A step of correcting the long-read assembly result, including using the library of the comparison result to correct the library of the split long-read assembly result to obtain a correction result;
[0023] A step of merging the correction results, including merging the correction results to obtain the corrected long-read assembly result.
[0024] It should be noted that the method of the present application for correcting long-read assembly results using short-read sequences first splits the long-read assembly results into smaller modules and then performs short-read Pilon correction. This not only significantly reduces the peak memory usage, but also can be executed in parallel when resources are sufficient, reducing the execution time. The error correction method of the present application has low requirements for device memory, does not require expensive large-memory nodes, nor does it need to lease expensive large-memory cloud servers. It not only improves the efficiency of correcting long-read assembly results using short-read sequences, but also reduces the error correction cost, laying a foundation for the further popularization and application of long-read sequencing.
[0025] It can be understood that one of the keys of the present application lies in pre-splitting the long-read assembly results. As for the specific splitting threshold, that is, the fixed threshold, it can be determined according to the long-read assembly results to be processed or the usage requirements. Other steps, such as library construction, alignment, error correction, etc., can refer to the prior art, except that the specific object of library construction in the present application is different. For example, in the library construction step of the long-read assembly results of the present application, the object of library construction is the split long-read assembly results, while the existing library construction 101 of the long-read assembly results, as Figure 2 shown, directly uses the long-read assembly results for library construction.
[0026] In one implementation of the present application, the fixed threshold is 100 Mb - 1000 Mb.
[0027] Preferably, the fixed threshold is 100 Mb.
[0028] It should be noted that the size of the genome is generally greater than 1 Gb and less than 100 Gb, and can be split according to the range of 100 Mb - 1000 Mb, and it is possible to complete genome error correction using a machine with a smaller peak memory usage. It can be understood that the fixed threshold can also be greater than 1000 Mb. Of course, the larger the threshold, the greater the corresponding peak memory required. It can be understood that the value of the fixed threshold actually depends on how many parts the long-read assembly results need to be split into. The fixed threshold is the size of each part. In principle, as long as the size of each part is limited to the extent that can be processed by the existing computer memory. Or, in order to further reduce memory occupancy, the fixed threshold can be designed to be smaller, which is not specifically limited here.
[0029] In an implementation of the present application, the long-read assembly results are segmented using a fixed threshold. Specifically, it includes setting a fixed threshold to segment the long-read assembly results. When the length of a contig reaches or exceeds the fixed threshold, it is regarded as one portion; when the length of a contig is less than the fixed threshold, it is combined with the adjacent contig in front or behind. When the cumulative length of the combined contigs reaches or exceeds the threshold, the combined contig is regarded as one portion; and so on. If the cumulative length of the last contig or the combined contig is less than the fixed threshold, it is also regarded as one portion.
[0030] It should be noted that the segmentation in the present application means that, based on a fixed threshold, segmentation is performed on the basis of one contig; that is, a single contig or multiple adjacent contigs that are just greater than or equal to the fixed threshold are regarded as one portion, and the remaining contigs, even if they are less than the fixed threshold size, are also regarded as one portion. For example, if a contig is greater than the fixed threshold, then this contig is directly regarded as one portion; if a contig is less than the fixed threshold, it needs to be combined with the adjacent contig. The number of contigs in the combination is such that the cumulative length after combination is just greater than or equal to the fixed threshold, and this combined contig is regarded as one portion; all contigs are segmented according to this principle; until the last contig. If the last contig is greater than or equal to the fixed threshold, it is directly regarded as one portion. If the last contig is less than the fixed threshold, it is also regarded as one portion. Therefore, in the long-read assembly result library after segmentation in the present application, except for the last portion which may be less than the fixed threshold, the remaining portions are all greater than or equal to the fixed threshold. For example, in an implementation of the present application, the 2.18 Gb long-read assembly results are segmented with a fixed threshold of 100 Mb, and 20 portions are obtained. The size of each portion is in sequence: 108 Mb, 114 Mb, 114 Mb, 110 Mb, 120 Mb, 109 Mb, 109 Mb, 109 Mb, 108 Mb, 108 Mb, 109 Mb, 108 Mb, 111 Mb, 109 Mb, 119 Mb, 108 Mb, 110 Mb, 108 Mb, 108 Mb, 46 Mb.
[0031] It should also be noted that when the present application performs segmentation using a fixed threshold, the contig is not truncated. The segmentation is based on the length of the contig or the cumulative length of the combined contigs reaching or exceeding the fixed threshold. Truncating the contig will affect the accuracy of the results.
[0032] Another aspect of the present application discloses an apparatus for correcting long-read assembly results using short-read sequences. The apparatus includes a cutting module, a library construction module for long-read assembly results, an alignment module, an alignment result matching module, an alignment result conversion and sorting module, a library construction module for alignment results, a long-read assembly result error correction module, and a merged error correction result module;
[0033] The cutting module includes means for splitting the long-read assembly result at a fixed threshold to obtain the split long-read assembly result;
[0034] The library construction module for long-read assembly results includes means for constructing a library for the split long-read assembly result to obtain a library of split long-read assembly results;
[0035] The alignment module includes means for aligning the short-read sequences to the library of split long-read assembly results to obtain an alignment result;
[0036] The alignment result matching module includes means for matching the alignment result obtained by the alignment module to the library of split long-read assembly results to obtain an alignment result containing matching information;
[0037] The alignment result conversion and sorting module includes means for converting the alignment result containing matching information into the bam format and sorting;
[0038] The library construction module for alignment results includes means for constructing a library for the sorted alignment result converted into the bam format to obtain a library of alignment results;
[0039] The long-read assembly result error correction module includes means for correcting the library of split long-read assembly results using the library of alignment results to obtain an error correction result;
[0040] The merged error correction result module includes means for merging the error correction results to obtain the error-corrected long-read assembly result.
[0041] It should be noted that the apparatus for correcting long-read assembly results using short-read sequences in the present application actually realizes each step in the method for correcting long-read assembly results using short-read sequences in the present application through each module; therefore, the specific limitations of each module can refer to the method for correcting long-read assembly results using short-read sequences in the present application. For example, the fixed threshold is 100 Mb, and specifically how to split, etc.
[0042] Another aspect of the present application discloses an apparatus for correcting long-read assembly results using short-read sequences. The apparatus includes a memory and a processor; wherein, the memory includes means for storing a program; the processor includes means for implementing the method for correcting long-read assembly results using short-read sequences in the present application by executing the program stored in the memory.
[0043] Another aspect of the present application discloses a computer-readable storage medium storing a program that can be executed by a processor to implement the method of using short read sequences to correct the long read assembly result in the present application.
[0044] Due to the above technical solutions, the beneficial effects of the present application are as follows:
[0045] The method and device for using short read sequences to correct the long read assembly result in the present application can use small memory nodes to correct the long read assembly result with short read sequences, greatly reducing the memory peak value, solving the problem that short read sequences cannot correct the long read assembly result without large memory nodes, and reducing the error correction cost of long assembly results; moreover, when resources are sufficient, the error correction method and device of the present application can also be executed in parallel, reducing the execution time and improving the error correction efficiency. Description of the Drawings
[0046] Figure 1 is the main process block diagram of the existing long read assembly;
[0047] Figure 2 is the process block diagram of the existing method of using short read sequences to correct the long read assembly result;
[0048] Figure 3 is the process block diagram of using short read sequences to correct the long read assembly result in the embodiment of the present application;
[0049] Figure 4 is a partial comparison result screenshot of the short read sequences aligned to the segmented long read assembly result library in the embodiment of the present application;
[0050] Figure 5 is an error screenshot of task interruption due to insufficient memory when using the existing method to correct the long read assembly result with short read sequences in the embodiment of the present application. Detailed Embodiments
[0051] The present application will be further described in detail below in conjunction with the drawings through specific embodiments. In the following embodiments, many detailed descriptions are provided to enable a better understanding of the present application. However, those skilled in the art can easily recognize that some of these features can be omitted in different situations, or can be replaced by other elements, materials, or methods. In some cases, some operations related to the present application are not shown or described in the specification to avoid overwhelming the core part of the present application with excessive descriptions. For those skilled in the art, it is not necessary to describe these related operations in detail, and they can fully understand the related operations based on the descriptions in the specification and general technical knowledge in the art.
[0052] Existing methods for correcting long-read assembly results using short-read sequences usually require large-memory nodes and renting computers with large memory (100 Gb), which greatly increases the cost of correcting long-read assembly results and is not conducive to the popularization and application of long-read sequencing.
[0053] The present application creatively proposes to split the long-read assembly result into smaller modules. When performing alignment and error correction, the short-read sequences sequentially perform alignment and error correction on each module, so that large-memory nodes are not required, and the short-read sequences can also correct the long-read assembly result.
[0054] According to the above inventive concept, the method for correcting long-read assembly results using short-read sequences in the present application, as Figure 3 shown, includes a cutting step 201, a library construction step 202 for the long-read assembly result, an alignment step 203, a matching step 204 for the alignment result, a conversion and sorting step 205 for the alignment result, a library construction step 206 for the alignment result, an error correction step 207 for the long-read assembly result, and a step 208 for merging the error correction results.
[0055] Among them, the cutting step 201 includes splitting the long-read assembly result at a fixed threshold to obtain the split long-read assembly result.
[0056] In one implementation of the present application, the long-read assembly result A is split at a fixed threshold to obtain multiple split long-read assembly results (A1, A2... An). Set a fixed threshold to split the long-read assembly result A. When the length of the last contig (contiguous sequence) and the cumulative length of the previous contigs reach or exceed the threshold, it is regarded as one portion. The cumulative length of the remaining contigs starts from 0. When the cumulative length of the last one reaches or exceeds the threshold, it is regarded as another portion, and so on. If the cumulative length of the last portion is less than the threshold, it is also defined as one portion An.
[0057] The library construction step 202 for the long-read assembly result includes constructing a library for the split long-read assembly result to obtain a library of the split long-read assembly result.
[0058] In one implementation of the present application, the command "bwaindex A" is used to construct a library for the long-read assembly result A in the cutting step 201.
[0059] The alignment step 203 includes aligning the short-read sequences to the library of the split long-read assembly result to obtain an alignment result.
[0060] In one implementation of the present application, short read sequences are aligned to the long read assembly result A to obtain the alignment result B. Specifically, the short read sequences NGS_1.fq and NGS_2.fq are aligned to the long read assembly result A using the command "bwa mem -MANGS_1.fq NGS_2.fq > B" to obtain the alignment result B.
[0061] The alignment result matching step 204 includes matching the alignment result obtained in the alignment step to the segmented long read assembly result library to obtain the alignment result containing the matching information.
[0062] In one implementation of the present application, the alignment result B matches the ContigID according to the segmented long read assembly results (A1, A2... An) to obtain the sub-alignment results (B1, B2... Bn). The method is as follows: (1) Use the keyword "bwa" as the splitting flag for the table header and content. The lines where the "bwa" splitting flag is located for different assembly results are different. For example Figure 4 is a partial screenshot of the alignment result B displayed using the command "less -SN B". The first 1 to 1317 lines are the table header of the alignment result B, and starting from the 1318th line is the alignment result content. The line numbers in the first column (1309 - 1323) are displayed by the "less -SN" command and do not actually exist. They are displayed specifically for the convenience of explaining this step. (2) The ContigID in the third column of the alignment result content is matched with the ContigID in the long read assembly results (A1, A2... An) to obtain the sub-alignment results (B1, B2... Bn), where the table header of B is retained in the sub-alignment results (B1, B2... Bn). For example Figure 4 the ContigID "000044F|arrow" in the third column of the 1318th line in [example] also exists in the Contig ID of the long read assembly result A1, so the result of the 1318th line is the alignment result content of B1, and so on.
[0063] The alignment result conversion and sorting step 205 includes converting the alignment result containing the matching information into the bam format and sorting it.
[0064] In an implementation of this application, the sub-comparison results are converted into the BAM format and sorted to obtain the sorted results (C1.bam, C2.bam... Cn.bam). The method is as follows: (1) For the sub-comparison results (B1, B2... Bn) obtained in the comparison result matching step 204, use the command "samtools view -bS Bi > Bi.bam" respectively to obtain the results after format conversion (B1.bam, B2.bam... Bn.bam), where i in the command Bi represents (1, 2... n). (2) Use the command "samtools sort Bi.bam Ci" to sort (B1.bam, B2.bam... Bn.bam) respectively to obtain the sorted results (C1.bam, C2.bam... Cn.bam), where i in the commands Bi and Ci represents (1, 2... n).
[0065] The library building step 206 for the comparison results includes building a library for the sorted comparison results converted into the BAM format to obtain a comparison result library.
[0066] In an implementation of this application, libraries are built for the sorted sub-comparison results (C1.bam, C2.bam... Cn.bam) to obtain the library building results (C1.bam.bai, C2.bam.bai... Cn.bam.bai). Use the command "samtools index Ci.bam" to build libraries for (C1.bam, C2.bam... Cn.bam) respectively, where i in the command Ci.bam represents (1, 2... n).
[0067] The error correction step 207 for the long read assembly results includes correcting the split long read assembly result library using the comparison result library to obtain the error correction results.
[0068] In an implementation of this application, the sub-long read assembly results (A1, A2... An) are corrected to obtain the corrected results (D1.fasta, D2.fasta... Dn.fasta). The command "java -Xmx50G -jar pilon-1.22.jar --genome Ai --output Di --outdir. / --diploid --fix snps,indels --threads num --flank 0 --frags Ci.bam" is used. Here, 50G in -Xmx50G is the preset memory parameter to be used, which can be adjusted according to the actual situation during execution; num is the number of threads. Generally, the larger it is set, the shorter the running time, but the required memory will be larger; pilon-1.22.jar is the program of this version, and it can be replaced if there is an updated version. Here, i in Ai, Di, and Ci.bam in the command represents (1, 2... n), Di is the prefix of the output file name, and "--outdir. / " means the results are output to the current directory, which can be changed if needed. The functions of "--diploid", "--fix snps,indels", "--flank 0" and more parameters can be viewed in the help instructions of the Pilon software.
[0069] Step 208 of merging the corrected results includes merging the corrected results to obtain the corrected long read assembly results.
[0070] In an implementation of this application, the corrected results (D1.fasta, D2.fasta... Dn.fasta) are merged to obtain the final corrected result A.pilon.fasta. The command "cat D1.fasta D2.fasta... Dn.fasta > A.pilon.fasta" is used. Here, the file names D1.fasta, D2.fasta... Dn.fasta are separated by a space "", and A.pilon.fasta is the file name of the final corrected result, which can be renamed if needed without affecting the final result.
[0071] Those skilled in the art can understand that all or part of the functions of the above method can be implemented in a hardware manner or in a computer program manner. When all or part of the functions in the above method are implemented in a computer program manner, the program can be stored in a computer-readable storage medium, and the storage medium can include: read-only memory, random access memory, magnetic disk, optical disk, hard disk, etc., and the above functions can be implemented by the computer executing the program. For example, the program is stored in the memory of the device, and when the processor executes the program in the memory, all or part of the above functions can be implemented. In addition, when all or part of the functions in the above embodiments are implemented in a computer program manner, the program can also be stored in storage media such as a server, another computer, magnetic disk, optical disk, flash drive or mobile hard disk, and saved to the memory of the local device by downloading or copying, or the system of the local device is updated in version. When the processor executes the program in the memory, all or part of the functions in the above method can be implemented.
[0072] Therefore, based on the method of the present application, the present application proposes a device for using short read sequences to correct the long read assembly results, and the device includes a cutting module, a long read assembly result library building module, a comparison module, a comparison result matching module, a comparison result conversion and sorting module, a comparison result library building module, a long read assembly result error correction module, and a merged error correction result module.
[0073] The cutting module includes means for splitting the long read assembly result at a fixed threshold to obtain the split long read assembly result.
[0074] The long read assembly result library building module includes means for building a library for the split long read assembly result to obtain a split long read assembly result library.
[0075] The comparison module includes means for aligning the short read sequences to the split long read assembly result library to obtain a comparison result.
[0076] The comparison result matching module includes means for matching the comparison result obtained by the comparison module to the split long read assembly result library to obtain a comparison result containing matching information.
[0077] The comparison result conversion and sorting module includes means for converting the comparison result containing matching information into the bam format and sorting.
[0078] The comparison result library building module includes means for building a library for the sorted comparison result converted into the bam format to obtain a comparison result library.
[0079] The long read assembly result error correction module includes means for correcting the split long read assembly result library using the comparison result library to obtain an error correction result.
[0080] A merged error correction result module, which is used to merge error correction results to obtain a corrected long read assembly result.
[0081] The device of the present application can implement the method of correcting the long read assembly result using short read sequences by the coordinated interaction of each module. In particular, each module of the device of the present application can implement the corresponding steps in the method of the present application, thereby realizing the automated correction of the long read assembly result using short read sequences.
[0082] In another implementation manner of the present application, a device for correcting a long read assembly result using short read sequences is also provided. The device includes a memory and a processor; the memory is used to store a program; the processor is used to implement the following method by executing the program stored in the memory: a cutting step, which includes cutting the long read assembly result with a fixed threshold to obtain a cut long read assembly result; a long read assembly result library building step, which includes building a library for the cut long read assembly result to obtain a cut long read assembly result library; a comparison step, which includes comparing the short read sequences with the cut long read assembly result library to obtain a comparison result; a comparison result matching step, which includes matching the comparison result obtained in the comparison step with the cut long read assembly result library to obtain a comparison result containing matching information; a comparison result conversion and sorting step, which includes converting the comparison result containing matching information into the bam format and sorting; a comparison result library building step, which includes building a library for the sorted comparison result converted into the bam format to obtain a comparison result library; a long read assembly result error correction step, which includes correcting the cut long read assembly result library using the comparison result library to obtain an error correction result; a merged error correction result step, which includes merging the error correction results to obtain a corrected long read assembly result.
[0083] In another implementation manner of the present application, a computer-readable storage medium is further provided. A program is stored in the storage medium, and the program can be executed by a processor to implement the following method: a cutting step, including cutting the long-read assembly result with a fixed threshold to obtain the cut long-read assembly result; a long-read assembly result library building step, including building a library for the cut long-read assembly result to obtain a cut long-read assembly result library; a comparison step, including aligning the short-read sequence to the cut long-read assembly result library to obtain a comparison result; a comparison result matching step, including matching the comparison result obtained in the comparison step to the cut long-read assembly result library to obtain a comparison result containing matching information; a comparison result conversion and sorting step, including converting the comparison result containing matching information into the bam format and sorting; a comparison result library building step, including building a library for the sorted comparison result converted into the bam format to obtain a comparison result library; a long-read assembly result error correction step, including correcting the cut long-read assembly result library with the comparison result library to obtain an error correction result; a combined error correction result step, including combining the error correction results to obtain a corrected long-read assembly result.
[0084] The present application will be further described in detail below through specific embodiments and drawings. The following embodiments are only for further illustrating the present application and should not be construed as limiting the present application.
[0085] Embodiment
[0086] In this example, the same long-read assembly result is corrected by using the existing conventional method and the improved method of this example respectively, and the advantages and effects of the improved method of this example compared with the conventional method are compared.
[0087] Among them, the existing conventional method, that is, the existing conventional method for correcting the long-read assembly result with short-read sequences, does not directly use the long-read assembly result for subsequent error correction without cutting. The improved method of this example is to cut the long-read assembly result and then use it for subsequent error correction.
[0088] In this example, the 199 Gb PacBio data in the Methods section of the article "Genome assembly of a tropical maize inbred line provides insights into structural variation and crop improvement" is assembled with the Falcon software and polished to obtain the assembly result A.fasta, with a size of about 2.18 Gb. This assembly result is used as the long-read assembly result for subsequent error correction processing.
[0089] The detailed link to the article is: www.nature.com / articles / s41588-019-0427-6.
[0090] Meanwhile, in this example, the following short-read data is downloaded:
[0091] sra-pub-src-2.s3.amazonaws.com / SRR8873351 / Clean_450_1.fq.gz
[0092] and sra-pub-src-2.s3.amazonaws.com / SRR8873351 / Clean_450_2.fq.gz
[0093] And for each of them, the first 59 Gb of data is intercepted to obtain 1.fq and 2.fq, that is, a total of 118 Gb of short-read data, which is used for error correction of the long-read assembly result A.fasta. The following uses Experiment 1 and Experiment 2 to illustrate the comparison effects of the implementation of the conventional method and the improved method in this example.
[0094] Experiment 1
[0095] As Figure 2 shown, in this experiment, the conventional method is used to perform pilon error correction on the aforementioned assembly result A.fasta to obtain the error correction result Pilon.fasta. The specific steps are as follows:
[0096] Step (1), building a library for the long assembly result 101, using the command "bwa index A.fasta" to build a library for the long-read assembly result to obtain the library building result.
[0097] Step (2), aligning the short-read sequences to the long assembly result 102, using the command "bwa mem -t80 -M A.fasta 1.fq 2.fq > B.sam" to align the short-read sequences 1.fq and 2.fq to the assembly result A.fasta to obtain the alignment result B.sam.
[0098] Step (3), converting the format of the comparison result 103, using the command "samtools view -bS B.sam > B.bam" to convert the B.sam format to B.bam to obtain the result B.bam.
[0099] Step (4), sorting the comparison result 104, using the command "samtools sort B.bam B" to sort the B.bam file and obtain the file with the same name B.bam.
[0100] Step (5), Build the library for the comparison result 105. Use the command "samtools index B.bam" to build the library and obtain the library building result B.bam.bai.
[0101] Step (6), Correct the errors in the long-length assembly result 106. Use the command "java -Xmx300G -jar pilon-1.22.jar --genome A.fasta --output Pilon --outdir. / --diploid --fix snps,indels --threads 100 --flank 0 --frags B.bam" to correct the errors in A.fasta and obtain the final error correction result Pilon.fasta. For this step task, in this experiment, we tried using 100G, that is, setting the parameter -Xmx100G, but reported the following Figure 5 error, and the task was interrupted because the memory was insufficient, that is, "java.lang.OutOfMemoryError". In addition, in this experiment, we also tried 200G, that is, setting the parameter -Xmx200G, but the task ran for more than 2 days without outputting results and without reporting errors. Later, after changing it to 300G in this experiment, the result was quickly obtained in 74 minutes.
[0102] The resources and time consumption of steps (1) to (6) of Experiment 1 are shown in Table 1.
[0103] Table 1 Resources and time consumption of each step of Experiment 1
[0104] Step Number of threads * Number of tasks Peak memory (Gb) Time (minutes) (1) 1*1 2.2 62 (2) 80*1 16 578 (3) 1*1 0.1 547 (4) 1*1 2 752 (5) 1*1 1 20 (6) 100*1 268 74
[0105] The results of Experiment 1 show that when using the existing conventional method of correcting the long-read assembly result with short-read sequences, a large memory node, such as 300G, is required to effectively correct the errors.
[0106] Experiment 2
[0107] In this experiment, the improved method of this example was used to perform pilon error correction on the same assembly result A.fasta of Experiment 1, and the error correction result Pilon.fasta was obtained. The steps are as follows:
[0108] Step (1), cutting step 201: Set a fixed threshold of 100MB to split the long-read assembly result A.fasta. When the length of the last contig and the cumulative length of the previous contigs reach or exceed this threshold of 100MB, it is regarded as one portion. The cumulative length of the remaining contigs starts from 0. When the cumulative length of the last contig reaches or exceeds the threshold of 100MB, it is regarded as another portion, and so on. If the cumulative length of the last portion is less than the threshold of 100MB, it is also defined as one portion. Finally, in this experiment, A.fasta is split into 20 portions, namely (A1.fasta, A2.fasta…A20.fasta), with sizes approximately (108Mb, 114Mb, 114Mb, 110Mb, 120Mb, 109Mb, 109Mb, 109Mb, 108Mb, 108Mb, 109Mb, 108Mb, 111Mb, 109Mb, 119Mb, 108Mb, 110Mb, 108Mb, 108Mb, 46Mb).
[0109] Step (2), library construction step 202 for the long-read assembly result: Use the command "bwa index A.fasta" to construct a library for the long-read assembly result. This step is the same as step (1) of Experiment 1. The difference is that the object of this experiment is the split long-read assembly result.
[0110] Step (3), alignment step 203: Use the command "bwa mem -t 80 -M A.fasta 1.fq 2.fq > B.sam" to align the short-read sequences 1.fq and 2.fq to the assembly result A.fasta to obtain the alignment result B.sam. This step is the same as step (2) of Experiment 1.
[0111] Step (4), alignment result matching step 204: The alignment result B.sam in step (3) matches the ContigID according to the split long-read assembly result (A1.fasta, A2.fasta…A20.fasta) in step (1) to obtain the sub-alignment results (B1.sam, B2.sam…B20.sam). The method is as follows: (1) Use the keyword "bwa" as the splitting flag for the header and content. For example Figure 4It is a partial screenshot of the comparison result B. Lines 1 to 1317 are the header of B.sam of the comparison result, and the comparison result content starts from line 1318. (2) For the third column ContigID in the comparison result content that matches the ContigID in the long read assembly results (A1.fasta, A2.fasta... A20.fasta) obtained by splitting in step (1), the sub-comparison results (B1.sam, B2.sam... B20.sam) are obtained, and the header of B.sam is retained in each of the sub-comparison results (B1.sam, B2.sam... B20.sam). For example Figure 4 The ContigID "000044F|arrow" in the third column of line 1318 in Figure 4 also exists in the Contig ID of the long read assembly result A8.fasta. Then the result of line 1318 is the comparison result content of B8.sam, and so on. Finally, the sub-comparison results (B1.sam, B2.sam... B20.sam) are obtained.
[0112] Step (5), comparison result conversion and sorting step 205, the sub-comparison results (B1.sam, B2.sam... B20.sam) are converted to the bam format and sorted to obtain the sorted results (C1.bam, C2.bam... C20.bam). The method is as follows: (1) Use the command "samtools view -bSBi>Bi.bam" for each of the sub-comparison results (B1.sam, B2.sam... B20.sam) to obtain the results after format conversion (B1.bam, B2.bam... B20.bam), where i in the command Bi represents (1, 2... 20). (2) Use the command "samtools sort Bi.bamCi" to sort each of (B1.bam, B2.bam... B20.bam) respectively to obtain the sorted results (C1.bam, C2.bam... C20.bam), where i in the commands Bi and Ci represents (1, 2... 20).
[0113] Step (6), comparison result library building step 206, use the command "samtools indexCi.bam" to build libraries for each of the sorted sub-comparison results (C1.bam, C2.bam... C20.bam) respectively to obtain the library building results (C1.bam.bai, C2.bam.bai... C20.bam.bai). Where i in the command Ci.bam represents (1, 2... 20).
[0114] Step (7), the long-read assembly result error correction step 207, use the command "java -Xmx30G -jar pilon-1.22.jar --genome Ai.fasta --output Di --outdir. / --diploid --fix snps,indels --threads 5 --flank 0 --frags Ci.bam" to correct the sub-long-read assembly results (A1.fasta, A2.fasta... A20.fasta), and obtain the error correction results (D1.fasta, D2.fasta... D20.fasta). Here, i in Ai, Di, and Ci.bam in the command represents (1, 2... 20), and the sub-error correction results (D1.fasta, D2.fasta... D20.fasta) are obtained.
[0115] Step (8), the merged error correction result step 208, use the command "cat D1.fasta D2.fasta... D20.fasta > A.pilon.fasta" to merge the error correction results (D1.fasta, D2.fasta... D20.fasta), and obtain the final error correction result A.pilon.fasta.
[0116] Table 2 Resource and time consumption of each step in Experiment 2
[0117] Step Number of threads * Number of tasks Peak memory (Gb) Average time (minutes) Longest time (minutes) (1) 1*1 3 2 2 (2) 1*1 2.2 62 62 (3) 80*1 16 578 578 (4) 1*20 2.2 564 614 (5) 1*20 1 120 132 (6) 1*20 1 2 3 (7) 5*20 26 13 15 (8) 1*1 1 1 1
[0118] After comparing and analyzing the results obtained by sorting the ContigID order of the final error correction result A.pilon.fasta in Experiment 2 and the final error correction result Pilon.fasta in Experiment 1, it is consistent with Pilon.fasta, which proves that the results obtained by the improved method in this example are reliable. In addition, according to the summary of Table 1 and Table 2, the resource consumption comparison between Experiment 1 and Experiment 2 is shown in Table 3. The total time consumption of Experiment 2 is 1407 minutes, which is 69.2% of the total time consumption of Experiment 1 (2033 minutes), significantly shortening the error correction time. The maximum memory peak value of Experiment 2 is 26 Gb, which is 9.7% of the maximum memory peak value of Experiment 1 (268 Gb), greatly reducing the memory peak value. This enables the improved method in this example to run on general small-memory nodes, such as computers with 32 Gb and 64 Gb of memory, and can meet the computer configuration status of most companies and universities.
[0119] Table 3 Resource consumption comparison between Experiment 1 and Experiment 2
[0120]
[0121]
[0122] In summary, the improved method in this example can use small memory nodes to correct the long-read assembly result Pilon with short-read sequences, significantly reducing the peak memory usage, solving the problem that many universities or companies cannot correct the long-read assembly result Pilon with short-read sequences due to the lack of large memory nodes, and also eliminating the need to lease cloud servers at high costs. Moreover, the improved method in this example can be executed in parallel when resources are sufficient, reducing the execution time and improving the error correction efficiency.
[0123] The above content is a further detailed description of the present application in combination with specific implementation manners, and it cannot be determined that the specific implementation of the present application is only limited to these descriptions. For those of ordinary skill in the technical field to which the present application belongs, without departing from the concept of the present application, several simple deductions or substitutions can still be made.
Claims
1. A method for correcting the results of long-read assembly using short-read sequences, characterized in that: comprising the following steps, a cutting step, including splitting the long-read assembly result at a fixed threshold to obtain the split long-read assembly result; a long-read assembly result library construction step, including constructing a library for the split long-read assembly result to obtain the split long-read assembly result library; an alignment step, including aligning the short-read sequences to the split long-read assembly result library to obtain an alignment result; an alignment result matching step, including matching the alignment result obtained in the alignment step to the split long-read assembly result library to obtain an alignment result containing matching information; an alignment result conversion and sorting step, including converting the alignment result containing matching information into the bam format and sorting it; an alignment result library construction step, including constructing a library for the sorted alignment result converted into the bam format to obtain an alignment result library; a long-read assembly result error correction step, including correcting the split long-read assembly result library using the alignment result library to obtain an error correction result; a step of merging the error correction results, including merging the error correction results to obtain the corrected long-read assembly result.
2. The method according to claim 1, characterized in that: The fixed threshold is 100 Mb - 1000 Mb.
3. The method according to claim 2, characterized in that: The fixed threshold is 100 Mb.
4. The method according to claim 1, characterized in that: The splitting of the long-read assembly result at a fixed threshold specifically includes setting a fixed threshold to split the long-read assembly result. When the length of a contig reaches or exceeds the fixed threshold, it is regarded as one portion; when the length of a contig is less than the fixed threshold, it is combined with the adjacent contig in front or behind. When the cumulative length of the combined contig reaches or exceeds the threshold, the combined contig is regarded as one portion; and so on. If the cumulative length of the last contig or combined contig is less than the fixed threshold, it is also regarded as one portion.
5. A device for correcting the results of long-read assembly using short-read sequences, characterized in that: including a cutting module, a long-read assembly result library construction module, an alignment module, an alignment result matching module, an alignment result conversion and sorting module, an alignment result library construction module, a long-read assembly result error correction module, and a step of merging the error correction results module; The cutting module includes means for splitting the long-read assembly result at a fixed threshold to obtain the split long-read assembly result; The long-read assembly result library construction module includes means for constructing a library for the split long-read assembly result to obtain the split long-read assembly result library; The alignment module includes means for aligning the short-read sequences to the split long-read assembly result library to obtain an alignment result; The alignment result matching module includes means for matching the alignment result obtained by the alignment module to the split long-read assembly result library to obtain an alignment result containing matching information; The alignment result conversion and sorting module includes means for converting the alignment result containing matching information into the bam format and sorting it; The alignment result library construction module includes means for constructing a library for the sorted alignment result converted into the bam format to obtain an alignment result library; The long-read assembly result error correction module includes correcting the segmented long-read assembly result library by using the comparison result library to obtain an error correction result; The merged error correction result module includes merging the error correction results to obtain a corrected long-read assembly result.
6. The device according to claim 5, characterized in that: The fixed threshold is 100MB.
7. The device according to claim 5, characterized in that: The segmentation of the long-read assembly result by using the fixed threshold specifically includes setting a fixed threshold to segment the long-read assembly result. When the length of a contig reaches or exceeds the fixed threshold, it is regarded as one portion; when the length of a contig is less than the fixed threshold, it is combined with the adjacent contig in front of or behind it. When the cumulative length of the combined contig reaches or exceeds the threshold, the combined contig is regarded as one portion; and so on. If the cumulative length of the last contig or the combined contig is less than the fixed threshold, it is also regarded as one portion.
8. A device for correcting the results of long-read assembly using short-read sequences, characterized in that: The device includes a memory and a processor; The memory includes storing a program; The processor includes implementing the method according to any one of claims 1-4 by executing the program stored in the memory.
9. A computer-readable storage medium, characterized in that: A program is stored in the storage medium, and the program can be executed by the processor to implement the method according to any one of claims 1-4.
Citation Information
Patent Citations
Method and system for assembling genomic sequence
CN104017883A
Splicing method and system of second generation and third generation genomic sequencing data combination
CN104951672A