Genome assembly method, system and equipment based on perl language and storage medium
Through a genome assembly method based on the Perl language, high-precision short reads collection and Hi-C sequencing data were used to correct the interaction intensity matrix of long reads, which solved the assembly problems of highly repetitive and large genomes and achieved high-quality genome assembly and continuity.
Patent Information
- Application Number
- CN202510720784.3
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-05-30
- Publication Date
- 2025-10-10
AI Technical Summary
In existing technologies, sequencing methods for highly repetitive and large genomes are economically expensive and suboptimal, making it difficult to achieve effective genome assembly.
A genome assembly method based on the Perl language was used. Through the segmentation of high-precision short reads sets, sequencing depth analysis and normal distribution model, combined with Hi-C sequencing data, the interaction intensity matrix of high-precision long reads was corrected, clustering and assembly were performed, and the influence of chromosomal repetitive sequences was reduced.
It improves the quality and continuity of genome assembly, reduces sequencing costs, and is suitable for the assembly needs of highly repetitive and large genomes.
Smart Images
Figure CN120766768A_ABST
Abstract
Description
Technical Field
[0001] The present application relates to the field of bioinformatics analysis technology, and in particular to a genome assembly method, system, device and storage medium based on the Perl language. Background Art
[0002] The genome contains all the genetic sequence information of an organism. With the advancement of sequencing technology, genome assembly techniques have become increasingly sophisticated, and more and more genomes have been deciphered. Currently, genome assembly methods often use third-generation sequencing (NGS) long reads to assemble a preliminary draft genome. Based on this draft, the interaction information between different chromosome segments is obtained through experimental and sequencing techniques, taking advantage of the fact that chromosomes generally have territoriality (i.e., sequences within the same chromosome are physically closer in space). By back-aligning the positions of sequencing reads on the draft genome, the interaction strengths between different contigs in the draft genome can be determined. Based on the decay pattern of interaction strengths along different sequence directions and the overall interaction strength between contigs, the spatial positions of two contigs on the chromosome can be inferred, and the genome can then be assembled at the chromosome level. Due to the increasing read lengths of NGS data, it is now possible to construct a complete genome map for some smaller genomes, that is, a telomere-to-telomere chromosome assembly. These ultra-long reads can resolve the complexity of regions where repetitive sequences are located in the genome, but because the economic cost of this sequencing method is still relatively high and the data generated by the sequencing machine is limited, it is still not an ideal method for genomes with high repetitiveness and large genomes. Summary of the Invention
[0003] In view of this, the present application provides a genome assembly method, system, device and storage medium based on the Perl language, aiming to improve the problem in the prior art that there is no ideal sequencing method for genomes with high repetitiveness and large genomes.
[0004] The embodiment of the present application is implemented as follows: a genome assembly method based on the Perl language comprises the following steps: Obtain a high-precision short Reads set, and split the 1 to SRL sites of each Read in the high-precision short Reads set into SRL-L+1 small fragments of length L, which are recorded as S set, where the elements are S i ; By S i The distribution of depth determines the average sequencing depth AD of the high-precision short read set, and the sequencing depth is calculated by the interval from AD to the inflection point that does not conform to the Poisson distribution data, which conforms to the standard deviation STD of the Poisson distribution data; A normal distribution model was used to fit the theoretical expected value of sequencing depth based on AD and STD. The upper limit Emin and lower limit Emax of the theoretical expected sequencing depth were defined based on the coverage area reaching 95% on both sides of the normal distribution mean. Based on the size of Emin and Emax, Si in the range from the upper limit Emin to the lower limit Emax is filtered out from the S set, which is recorded as the filtered Sf set; Get the first effective Reff i ; Get a high-precision long Reads set and record each Read in the high-precision long Reads set as HR i , based on the first effective Reff i Get the second effective Reffi and determine the second effective Reff i The corresponding high-precision long Reads in the genome, if the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN; According to the level of Reff_set_filterN, keep the same as HR i Corresponding to the high-precision long Reads with the largest number of Sfs, the high-precision long Reads corresponding to Reff_set_filter are recorded as Reff_set_filter_best; According to Reff_set_filter_best, get any two HRs i The original interaction strength two-dimensional matrix ForceMatrix; Use any value v in ForceMatrix ij The lengths of the two corresponding high-precision long reads are corrected to obtain the corrected high-precision long read interaction strength two-dimensional matrix NormalForceMatrix; Split the two-dimensional matrix NormalForceMatrix to obtain high-precision clustering results of long reads; According to the clustering results of high-precision long reads, assembly is performed to obtain the final results.
[0005] Optionally, in some embodiments of the present application, a Hi-C sequencing Reads set is obtained, and the 1 to HRL sites of each Read in the Hi-C sequencing Reads set are divided into HRL-L+1 subsets to obtain the Sf of each Read. iSet, the number of sets is NSf i , Sf i The number of occurrences in the Sf set is recorded as ENSf i , use ENSf i Divide by NSf i , and obtain the ratio Re. When Re is greater than or equal to 0.8, it is considered that the read in the Hi-C sequencing reads set is suitable for identifying the clustering relationship between different reads in the high-precision long reads set, and is recorded as the first effective Reff i and / or Get a high-precision long Reads set and record each Read in the high-precision long Reads set as HR i , HR i 1 to LRL sites are divided into LRL-L+1 small segments of length L and recorded as set HRSET i and / or According to the different HRs corresponding to the forward and reverse reads in Reff_set_filter_best i Count the cases, record them as Hai, Hbi, and iterate them pair by pair to get any two HRs. i The two-dimensional matrix ForceMatrix of the raw interaction strength between them; and / or The two high-precision long reads are recorded as Hifi1 and Hifi2; the lengths of the two high-precision long reads are recorded as LengthHifi1 and LengthHifi2, and the correction formula is v ij *1e8 / (LengthHifi1+LengthHifi2), according to the correction formula, the corrected high-precision long reads interaction strength two-dimensional matrix NormalForceMatrix is obtained.
[0006] Optionally, in some embodiments of the present application, a high-precision short reads set of the target species is obtained by high-throughput sequencing; and / or Obtaining a Hi-C sequencing read set of the target species through high-throughput sequencing; and / or High-throughput sequencing is used to obtain a high-precision long read collection of the target species.
[0007] Optionally, in some embodiments of the present application, the first valid Reff is recorded by a hash data structure Hash-Reff i The number of elements in the Sf set ENSf i , ENSf iSort from large to small, and the top 95% of the values are recorded as the second effective Reffi. The second effective Reffi is calculated in turn with the intersection of HRSTi to obtain the number of elements Ninter contained in the intersection. If the number Ninter is less than half of the second Reffi length HRL, it is recorded as 0, and the rest are counted according to the actual number. If the high-precision long reads corresponding to the forward and reverse Hi-C sequencing reads come from the same line, then the second effective Reffi is considered i It does not provide effective information for clustering high-precision long Reads; if not, the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN.
[0008] Optionally, in some embodiments of the present application, during the pair-by-pair iteration process, if HR i If it is greater than the existing Hai, Hbi, the different HRs corresponding to the forward and reverse reads will be i replace.
[0009] Optionally, in some embodiments of the present application, the two-dimensional matrix NormalForceMatrix is split into a list of corrected intensities between any two high-precision long reads, the first and second columns represent two high-precision long reads a and b, and the third column represents the interaction intensity v; the clustering process first checks v to obtain a descending order of v. When the v value is the largest, the corresponding two high-precision long reads are first picked out to create a cluster and recorded, and then v is checked in sequence. If the Reads ID in the row has been recorded, it is included in the clustered group. If it has not been recorded, a new cluster is created and recorded again to finally obtain the clustering result of the high-precision long reads.
[0010] Optionally, in some embodiments of the present application, based on the clustering results of high-precision long reads, a certain type of high-precision long reads sequencing data in the clustering results is used as input, and the sequencing data of this type of clustering results are assembled using HiFiasm software, and then the high-precision long reads of other clustering results are assembled using the same method, and finally the results after merging different clustering results are merged to obtain the final assembly result.
[0011] Accordingly, the present application also provides a genome assembly system based on the Perl language, including: The high-precision short reads collection module is used to obtain a high-precision short reads collection. The 1 to SRL sites of each read in the high-precision short reads collection are divided into SRL-L+1 small fragments of length L, which are recorded as S sets. The elements of the S i ; AD and STD modules for S i The distribution of depth determines the average sequencing depth AD of the high-precision short read set, and the sequencing depth is calculated by the interval from AD to the inflection point that does not conform to the Poisson distribution data, which conforms to the standard deviation STD of the Poisson distribution data; The upper limit Emin and lower limit Emax position determination module is used to fit the theoretical expected value of sequencing depth based on AD and STD using a normal distribution model, and define the upper limit Emin and lower limit Emax of the theoretical expectation of sequencing depth based on the region where the coverage area on both sides of the normal distribution mean reaches 95%; The module for obtaining the filtered Sf set is used to filter out Si within the range of the upper limit Emin to the lower limit Emax from the S set based on the size of Emin and Emax, and record it as the filtered Sf set; First effective Reff i Get module, used to get the first effective Reff i ; Reff_set_filterN is a module used to obtain a high-precision long Reads set, and each Read in the high-precision long Reads set is recorded as HR i , based on the first effective Reff i Get the second effective Reffi and determine the second effective Reff i The corresponding high-precision long Reads in the genome, if the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN; Reads is the module for obtaining Reff_set_filter_best, which is used to retain the value of HR according to the level of Reff_set_filterN. i Corresponding to the high-precision long Reads with the largest number of Sfs, the high-precision long Reads corresponding to Reff_set_filter are recorded as Reff_set_filter_best; ForceMatrix acquisition module, used to obtain any two HRs according to Reff_set_filter_best ia raw interaction strength two-dimensional matrix ForceMatrix between the two high-precision long Reads; a NormalForceMatrix obtaining module configured to obtain any value v in the ForceMatrix ij correct the lengths of the two high-precision long Reads to obtain a corrected high-precision long Reads interaction strength two-dimensional matrix NormalForceMatrix; a high-precision long Reads clustering result obtaining module configured to split the two-dimensional matrix NormalForceMatrix to obtain a high-precision long Reads clustering result; a final result obtaining module configured to assemble a final result according to the high-precision long Reads clustering result.
[0012] Correspondingly, the embodiments of the present application further provide a computer device, comprising a storage and a processor, the storage stores a computer program, and the computer program is executed by the processor to make the processor execute the steps of the above method.
[0013] Correspondingly, the embodiments of the present application further provide a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to make the processor execute the steps of the above method.
[0014] In the present application, AD is Average depth, average sequencing depth; STD is standard error, standard deviation.
[0015] In the present application, a series of bioinformatics methods and mathematical algorithms are adopted, and the sequences measured by a single chromosome are classified into a group. Such processing can reduce the impact of other chromosome repetitive sequences (mainly from telomeres and centromeres) on assembly, thereby reducing the complexity of the assembly process of each chromosome, effectively increasing the sequence and continuity of a single chromosome (reflected in the average length of the assembled sequence), deeply utilizing the effective information of the data obtained by the sequencing technology, and then obtaining a better genome map in assembly quality. BRIEF DESCRIPTION OF DRAWINGS
[0016] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.
[0017] Figure 1 is a flowchart of the genome assembly method based on perl language in the present application. DETAILED DESCRIPTION
[0018] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of this application. In addition, it should be understood that the specific embodiments described herein are only used to illustrate and explain the present application and are not used to limit the present application.
[0019] In this application, "and / or" describes the association relationship between associated objects, indicating that three relationships can exist. For example, A and / or B can mean: A exists alone, A and B exist at the same time, and B exists alone. A and B can be singular or plural.
[0020] In this application, "at least one" means one or more, and "plurality" means two or more. "One or more," "at least one of the following," or similar expressions refer to any combination of these items, including any combination of single or plural items. For example, "at least one of a, b, or c" or "at least one of a, b and c" can all mean: a, b, c, ab (i.e., a and b), ac, bc, or abc, where a, b, and c can each be single or plural.
[0021] Various embodiments of the present application may be presented in the form of a range; it should be understood that the description in the form of a range is only for convenience and brevity and should not be understood as a hard limitation on the scope of the present application; therefore, the range description should be considered to have specifically disclosed all possible sub-ranges and single numbers within the range. For example, the description of a range from 1 to 6 should be considered to have specifically disclosed sub-ranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6, etc., as well as single numbers within the range, such as 1, 2, 3, 4, 5 and 6, which applies regardless of the range. In addition, whenever a numerical range is indicated herein, it is meant to include any cited number (fractional or integer) within the indicated range.
[0022] The technical solution of this application is as follows: See also Figure 1 The present invention provides a genome assembly method based on the Perl language, comprising the following steps: S01. Obtain a high-precision short reads set. Split the 1 to SRL sites of each read in the high-precision short reads set into SRL-L+1 small segments of length L, which are recorded as S set. The elements of the S set are S i , i is the number of times put in; S02, through S i The distribution of depth determines the average sequencing depth AD of the high-precision short read set, and the sequencing depth is calculated by the interval from AD to the inflection point that does not conform to the Poisson distribution data, which conforms to the standard deviation STD of the Poisson distribution data; S03. Use the normal distribution model to fit the theoretical expected value of sequencing depth based on AD and STD, and define the upper limit Emin and lower limit Emax of the theoretical expected sequencing depth based on the coverage area reaching 95% on both sides of the normal distribution mean; S04. Based on the size of Emin and Emax, filter out Si in the range from the upper limit Emin to the lower limit Emax from the S set, and record it as the filtered Sf set; S05. Get the first effective Reff i ; S06. Obtain a high-precision long Reads set and record each Read in the high-precision long Reads set as HR i , based on the first effective Reff i Get the second effective Reffi and determine the second effective Reff i The corresponding high-precision long Reads in the genome, if the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN; S07, according to the level of Reff_set_filterN, retain the same as HR i Corresponding to the high-precision long Reads with the largest number of Sfs, the high-precision long Reads corresponding to Reff_set_filter are recorded as Reff_set_filter_best; S08. According to Reff_set_filter_best, obtain any two HRs i The original interaction strength two-dimensional matrix ForceMatrix; S09, use any value v in ForceMatrix ijThe length of the two high-precision long reads is corrected to obtain a corrected high-precision long read interaction strength two-dimensional matrix NormalForceMatrix. S10, the two-dimensional matrix NormalForceMatrix is split to obtain a clustering result of the high-precision long reads; S11, according to the clustering result of the high-precision long reads, assembling to obtain a final result.
[0023] In the S01, In some embodiments, the target species high-precision short read set is obtained by high-throughput sequencing.
[0024] In some embodiments, 17≤L≤63, and L is an odd number.
[0025] It can be understood that the high-precision (error rate of each site is less than 1%) short read set (sequencing length is SRL) refers to that the read length of the genome second-generation sequencing is generally SRL=150bp, and the principle is to randomly take a sequence with a read length of SRL in the genome, but due to the existence of a certain error rate NGSe in sequencing, a certain error rate will exist in the sequencing result. The economic cost of second-generation sequencing is low, and the second-generation sequencing of the genome is required to reach at least 30 times of depth. Therefore, the expected depth of each base site in the genome is Ed=30; since the genome sequence is composed of ATCG four bases, in order to reflect the difference of each base type of each site, 1 to SRL sites in each read are divided into SRL-L+1 small fragments of L length, and the elements in the set S are S i .
[0026] Exemplarily, the sequencing read sequence of SRL=8bp is ATCGGCTA, which can be divided into four small fragments according to the length of L=5, which are S1=ATCGG, S2=TCGGC, S3=CGGCT, and S4=GGCTA. The set S contains S1 to S4 elements.
[0027] In the S02, It can be understood that the size of the general genome can reach 10e7~10e8 bp in length, or even longer, so that the set S contains a large number of non-repeated elements, and the number of 10e7~10e8 can also be obtained by sampling S i The probability P of sampling i is very small, so it can be considered that all S i depth conforms to the Poisson distribution model, but due to the existence of a certain error rate in sequencing and the existence of repeated regions in the genome and other biological characteristics, there must be a part of Si By S i The depth distribution determines the average sequencing depth AD of the second-generation sequencing data, and the standard deviation STD of the data that conforms to the Poisson distribution is determined by the distance from AD to the inflection point of the data that does not conform to the Poisson distribution.
[0028] In said S03: It can be understood that the normal distribution model is the asymptotic version of the Poisson distribution. After determining AD and STD, the normal distribution model is used to fit the theoretical expected value of the sequencing depth again.
[0029] In said S04: It is understandable that Sf contains the expected S i , these S i Minimally affected by genomic repetitive sequences and sequencing errors.
[0030] In said S05: In some embodiments, a Hi-C sequencing reads set is obtained, and the 1 to HRL sites of each read in the Hi-C sequencing reads set are divided into HRL-L+1 subsets to obtain the Sf of each read. i Set, the number of sets is NSf i , Sf i The number of occurrences in the Sf set is recorded as ENSf i , use ENSf i Divide by NSf i , and obtain the ratio Re. When Re is greater than or equal to 0.8, it is considered that the read in the Hi-C sequencing reads set is suitable for identifying the clustering relationship between different reads in the high-precision long reads set, and is recorded as the first effective Reff i .
[0031] Furthermore, a Hi-C sequencing reads set of the target species was obtained by high-throughput sequencing.
[0032] It can be understood that when obtaining a Hi-C sequencing reads set, the sequencing length is HRL.
[0033] Furthermore, the 1 to LRL sites of each read in the Hi-C sequencing reads set were mapped using the perl program.
[0034] Furthermore, a high-precision long reads collection of the target species was obtained through high-throughput sequencing.
[0035] It can be understood that in the high-precision long reads collection, the error rate of each site is less than 1%.
[0036] It can be understood that when a high-precision long read set is obtained, the sequencing length is LRL.
[0037] In said S06: In some embodiments, a high-precision long Reads set is obtained, and each Read in the high-precision long Reads set is recorded as HR i , HR i 1 to LRL sites are divided into LRL-L+1 small segments of length L and recorded as set HRSET i .
[0038] Furthermore, the first valid Reff is recorded by the hash data structure Hash-Reff i The number of elements in the Sf set ENSf i , ENSf i Sort from large to small, and the top 95% of the values are recorded as the second effective Reffi. The second effective Reffi is calculated in turn with the intersection of HRSTi to obtain the number of elements Ninter contained in the intersection. If the number Ninter is less than half of the second Reffi length HRL, it is recorded as 0, and the rest are counted according to the actual number. If the high-precision long reads corresponding to the forward and reverse Hi-C sequencing reads come from the same line, then the second effective Reffi is considered i It does not provide effective information for clustering high-precision long Reads; if not, the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN.
[0039] Furthermore, a high-precision long reads collection of the target species was obtained through high-throughput sequencing.
[0040] It is understandable that if N inter The number is lower than Reff i Half of the length HRL is recorded as 0, and the rest is counted according to the actual number to determine Reff i The corresponding high-precision long reads in the genome.
[0041] In said S08: In some embodiments, the different HRs corresponding to the forward and reverse reads in Reff_set_filter_best are used. i Count the cases, record them as Hai, Hbi, and iterate them pair by pair to get any two HRs. iThe original interaction strength two-dimensional matrix ForceMatrix between them.
[0042] Furthermore, in the pairwise iteration process, if HR i If it is greater than the existing Hai, Hbi, the different HRs corresponding to the forward and reverse reads will be i replace.
[0043] For example, L1 and R1, L2 and R2 have the following correspondence: L1 corresponds to high-precision long Reads 1, R1 corresponds to high-precision long Reads 2, L2 corresponds to high-precision long Reads 1, R2 corresponds to high-precision long Reads 2, then the interaction strength Fab between high-precision long Reads 1 and high-precision long Reads 2 is recorded as 2.
[0044] In some embodiments, in the two-dimensional matrix ForceMatrix, both the horizontal and vertical directions represent the IDs of all high-precision long reads, and the diagonal elements in the matrix correspond to the interaction strength of the same high-precision long reads; It can be understood that the Hi-C sequencing reads set represents a sequencing method that carries the interaction strength between the slices in different chromosomes, and the high-precision (error rate of each site is less than 1%) long reads set generally refers to HiFi sequencing.
[0045] It can be understood that in the two-dimensional matrix ForceMatrix, both the horizontal and vertical directions represent the IDs of all high-precision long reads. At this time, each value in the matrix can correspond to any two high-precision long reads, and the diagonal elements in the matrix correspond to the interaction strength of the same high-precision long reads. Since the Hi-C sequencing reads that fall into the same high-precision long read have been removed in step S06, they are all 0.
[0046] In said S09: It is understandable that since the lengths of reads in HiFi sequencing data are not consistent, the two-dimensional interaction strength matrix ForceMatrix has a certain bias. Specifically, the longer the two high-precision long reads, the larger the ForceMatrix value may be, but their natural physical positions are not close.
[0047] In some embodiments, two high-precision long reads are recorded as Hifi1 and Hifi2; the lengths of the two high-precision long reads are recorded as LengthHifi1 and LengthHifi2, and the correction formula is v ij*1e8 / (LengthHifi1+LengthHifi2), according to the correction formula, the corrected high-precision long reads interaction strength two-dimensional matrix NormalForceMatrix is obtained.
[0048] In the S10: In some embodiments, the two-dimensional matrix NormalForceMatrix is split into a list of corrected intensities between any two high-precision long reads, where the first and second columns represent two high-precision long reads a and b, and the third column represents the interaction intensity v; the clustering process first checks v to obtain a descending order of v. When the v value is the largest, the corresponding two high-precision long reads are first picked out to create a cluster and recorded, and then v is checked sequentially. If the Reads ID in the row has been recorded, it is included in the clustered group. If it has not been recorded, a new cluster is created and recorded again to finally obtain the clustering result of the high-precision long reads.
[0049] For example, the first line is 1 2 10, the second line is 1 4 11; the third line is 3 5 9; first sort by v from large to small, the result after sorting is, for example, the first line is 1 4 11, the second line is 1 2 10; the third line is 3 5 9; first create cluster g1; and use hash Hash_u to record that reads 1 and 4 have been used. At this time, check the second line; there are 3 situations in the second line. If the first and second lines are 2 3 10, it means that two reads a=2 and b=3 have not been recorded by Hash_u, so a new cluster g2 is created; if the second and second lines are 1 2 10, it means that two reads a=1 and b=2, at this time a is recorded by g, but b is not recorded by Hash_u, then b is classified into cluster g1 where a is; if the third and second lines are 2 410, both reads (a=2, b=4) are recorded, then all reads in the cluster where b is located are classified into the cluster where a is located, and finally the clustering of all HiFi Reads is completed.
[0050] Understandable, ij Where i and j are the i-th row and j-th column of the matrix.
[0051] In said S11: In some embodiments, based on the clustering results of high-precision long reads, a certain type of high-precision long reads sequencing data in the clustering results is used as input, and the sequencing data of this type of clustering results are assembled using HiFiasm software, and then the high-precision long reads of other clustering results are assembled using the same method, and finally the results of merging different clustering results are merged to obtain the final assembly result.
[0052] For example, g1: 1, 2, 3, g2: 4, 5, 6; based on the read IDs contained in g1 and g2, g1.fa and g2.fa are extracted from the total reads, representing the sequencing data of different clustering results, respectively. g1.fa and g2.fa are used as input, and HiFiasm software is used to assemble g1.fa and g2.fa separately. Then, g1.fa and g2.fa are merged into gcombine.fa as the final result of the divide-and-conquer assembly method for different chromosomes.
[0053] In a second aspect, the present invention provides a genome assembly system based on the Perl language, comprising: The high-precision short reads collection module is used to obtain a high-precision short reads collection. The 1 to SRL sites of each read in the high-precision short reads collection are divided into SRL-L+1 small fragments of length L, which are recorded as S sets. The elements of the S i ; AD and STD modules for S i The distribution of depth determines the average sequencing depth AD of the high-precision short read set, and the sequencing depth is calculated by the interval from AD to the inflection point that does not conform to the Poisson distribution data, which conforms to the standard deviation STD of the Poisson distribution data; The upper limit Emin and lower limit Emax position determination module is used to fit the theoretical expected value of sequencing depth based on AD and STD using a normal distribution model, and define the upper limit Emin and lower limit Emax of the theoretical expectation of sequencing depth based on the region where the coverage area on both sides of the normal distribution mean reaches 95%; The module for obtaining the filtered Sf set is used to filter out Si within the range of the upper limit Emin to the lower limit Emax from the S set based on the size of Emin and Emax, and record it as the filtered Sf set; First effective Reff i Get module, used to get the first effective Reff i ; Reff_set_filterN is a module used to obtain a high-precision long Reads set, and each Read in the high-precision long Reads set is recorded as HR i , based on the first effective Reff i Get the second effective Reffi and determine the second effective Reff i The corresponding high-precision long Reads in the genome, if the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reffi The number of overlapping Sf between different high-precision long reads Reff_set_filterN; Reads is the module for obtaining Reff_set_filter_best, which is used to retain the value of HR according to the level of Reff_set_filterN. i Corresponding to the high-precision long Reads with the largest number of Sfs, the high-precision long Reads corresponding to Reff_set_filter are recorded as Reff_set_filter_best; ForceMatrix acquisition module, used to obtain any two HRs according to Reff_set_filter_best i The original interaction strength two-dimensional matrix ForceMatrix; NormalForceMatrix obtains the module, using any value v in ForceMatrix ij The lengths of the two corresponding high-precision long reads are corrected to obtain the corrected high-precision long read interaction strength two-dimensional matrix NormalForceMatrix; The module for obtaining the clustering results of high-precision long reads is used to split the two-dimensional matrix NormalForceMatrix to obtain the clustering results of high-precision long reads; The final result acquisition module is used to assemble the final results based on the clustering results of high-precision long reads.
[0054] In the high-precision short reads collection module: In some embodiments, a high-precision short read collection of a target species is obtained by high-throughput sequencing.
[0055] In some embodiments, 17≤L≤63, and L is an odd number.
[0056] It can be understood that the high-precision (error rate of each site is less than 1%) short Reads set (sequencing length is SRL) refers to the second-generation sequencing of the genome, which is generally SRL=150bp in length. The principle is to randomly take a sequence with a read length of SRL in the genome, but due to a certain error rate in sequencing NGSe, it will lead to a certain error rate in the sequencing results. The economic cost of second-generation sequencing is low, and this step requires the second-generation sequencing of the genome to reach a depth of at least 30 times. Therefore, the expected depth of each base site in the genome is Ed=30; since the genome sequence is composed of four bases of ATCG, in order to reflect the difference of each base type at each site, the 1 to SRL sites in each Read are divided into SRL-L+1 small fragments of L length, recorded as S set, the elements of which are S i .
[0057] For example, a sequence of ATCGGCTA with a sequencing read sequence of SRL=8bp can be divided into four small fragments according to the length of L=5, namely S1=ATCGG, S2=TCGGC, S3=CGGCT, and S4=GGCTA, and the S set contains S1 to S4 elements.
[0058] In the AD and STD modules: It is understandable that the size of a general genome can reach 10e7~10e8 bp in length, or even longer. Therefore, if the set S contains a large number of non-repeated elements, 10e7~10e8 numbers can be obtained. i The probability of winning P i is very small, so it can be assumed that all S i The depth of S is consistent with the Poisson distribution model, but due to the certain error rate in sequencing and the biological characteristics of the genome such as repeated regions, there must be a part of S that does not conform to the Poisson distribution model. i By S i The depth distribution determines the average sequencing depth AD of the second-generation sequencing data, and the standard deviation STD of the data that conforms to the Poisson distribution is determined by the distance from AD to the inflection point of the data that does not conform to the Poisson distribution.
[0059] In the upper limit Emin and lower limit Emax position determination module: It can be understood that the normal distribution model is the asymptotic version of the Poisson distribution. After determining AD and STD, the normal distribution model is used to fit the theoretical expected value of the sequencing depth again.
[0060] In the module for obtaining the filtered Sf set: It is understandable that Sf contains the expected S i , these S iThe impact of genomic repeats and sequencing errors is extremely low.
[0061] The first effective Reff i In the obtaining module: In some embodiments, a set of Hi-C sequencing reads is obtained, and 1 to HRL positions of each read in the set of Hi-C sequencing reads are divided into HRL-L+1 subsets, to obtain Sf i sets, and the number of sets is NSf i The Sf i The number of occurrences of Sf in the set is recorded as ENSf i ENSf i is divided by NSf i , to obtain a ratio Re, and when Re is greater than or equal to 0.8, it is considered that the read in the set of Hi-C sequencing reads is suitable for identifying the clustering relationship between different reads in the set of high-precision long reads, and is recorded as the first effective Reff i .
[0062] Further, a set of high-precision long reads of the target species is obtained by high-throughput sequencing.
[0063] It can be understood that when the set of Hi-C sequencing reads is obtained, the sequencing length is HRL.
[0064] Further, 1 to LRL positions of each read in the set of Hi-C sequencing reads are divided by a perl program.
[0065] Further, a set of high-precision long reads of the target species is obtained by high-throughput sequencing.
[0066] It can be understood that in the set of high-precision long reads, the error rate of each position is less than 1%.
[0067] It can be understood that when the set of high-precision long reads is obtained, the sequencing length is LRL.
[0068] In the Reff_set_filterN obtaining module: In some embodiments, a set of high-precision long reads is obtained, and each read in the set of high-precision long reads is recorded as HR i 1 to LRL positions of HR i are divided into LRL-L+1 small fragments of L length, recorded as a set HRSET i .
[0069] Further, the first effective Reff is recorded by a hash data structure Hash-Reffi The number of elements in the Sf set ENSf i , ENSf i Sort from large to small, and the top 95% of the values are recorded as the second effective Reffi. The second effective Reffi is calculated in turn with the intersection of HRSTi to obtain the number of elements Ninter contained in the intersection. If the number Ninter is less than half of the second Reffi length HRL, it is recorded as 0, and the rest are counted according to the actual number. If the high-precision long reads corresponding to the forward and reverse Hi-C sequencing reads come from the same line, then the second effective Reffi is considered i It does not provide effective information for clustering high-precision long Reads; if not, the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN.
[0070] Furthermore, a high-precision long reads collection of the target species was obtained through high-throughput sequencing.
[0071] It is understandable that if N inter The number is lower than Reff i Half of the length HRL is recorded as 0, and the rest is counted according to the actual number to determine Reff i The corresponding high-precision long reads in the genome.
[0072] The ForceMatrix acquisition module: In some embodiments, the different HRs corresponding to the forward and reverse reads in Reff_set_filter_best are used. i Count the cases, record them as Hai, Hbi, and iterate them pair by pair to get any two HRs. i The original interaction strength two-dimensional matrix ForceMatrix between them.
[0073] Furthermore, in the pairwise iteration process, if HR i If it is greater than the existing Hai, Hbi, the different HRs corresponding to the forward and reverse reads will be i replace.
[0074] Exemplarily, L1 and R1, L2 and R2 have the following corresponding relationship: L1 corresponds to high-precision long Reads 1, R1 corresponds to high-precision long Reads 2, L2 corresponds to high-precision long Reads 1, and R2 corresponds to high-precision long Reads 2. The interaction strength Fab of high-precision long Reads 1 and high-precision long Reads 2 is recorded as 2.
[0075] In some embodiments, in the two-dimensional matrix ForceMatrix, the horizontal and vertical directions both represent the IDs of all high-precision long Reads, and the diagonal elements in the matrix correspond to the interaction strength of the same high-precision long Reads. It can be understood that the Hi-C sequencing Reads set represents a sequencing method carrying the interaction strength between fragments in different chromosomes, and the high-precision (error rate of each site is less than 1%) long Reads set generally refers to HiFi sequencing.
[0076] It can be understood that in the two-dimensional matrix ForceMatrix, the horizontal and vertical directions both represent the IDs of all high-precision long Reads, and each value in the matrix can correspond to any two high-precision long Reads, and the diagonal elements in the matrix correspond to the interaction strength of the same high-precision long Reads. Since the Hi-C sequencing Reads falling into the same high-precision long Reads have been removed in step S06, they are all 0.
[0077] In the NormalForceMatrix obtaining module: It can be understood that since the lengths of the Reads in the HiFi sequencing data are not consistent, the interaction strength two-dimensional matrix ForceMatrix has certain preference, which is specifically manifested in that the longer the lengths of the two high-precision long Reads are, the larger the ForceMatrix value may be, but the actual physical positions thereof are not close.
[0078] In some embodiments, the two high-precision long Reads are recorded as Hifi1 and Hifi2, the lengths of the two high-precision long Reads are recorded as LengthHifi1 and LengthHifi2, and the correction formula is v ij *1e8 / ( LengthHifi1+LengthHifi2), and the corrected high-precision long Reads interaction strength two-dimensional matrix NormalForceMatrix is obtained according to the correction formula.
[0079] In the high-precision long Reads clustering result obtaining module: In some embodiments, the two-dimensional matrix NormalForceMatrix is split into a list of corrected intensities between any two high-precision long reads, where the first and second columns represent two high-precision long reads a and b, and the third column represents the interaction intensity v; the clustering process first checks v to obtain a descending order of v. When the v value is the largest, the corresponding two high-precision long reads are first picked out to create a cluster and recorded, and then v is checked sequentially. If the Reads ID in the row has been recorded, it is included in the clustered group. If it has not been recorded, a new cluster is created and recorded again to finally obtain the clustering result of the high-precision long reads.
[0080] For example, the first line is 1 2 10, the second line is 1 4 11; the third line is 3 5 9; first sort by v from large to small, the result after sorting is, for example, the first line is 1 4 11, the second line is 1 2 10; the third line is 3 5 9; first create cluster g1; and use hash Hash_u to record that reads 1 and 4 have been used. At this time, check the second line; there are 3 situations in the second line. If the first and second lines are 2 3 10, it means that two reads a=2 and b=3 have not been recorded by Hash_u, so a new cluster g2 is created; if the second and second lines are 1 2 10, it means that two reads a=1 and b=2, at this time a is recorded by g, but b is not recorded by Hash_u, then b is classified into cluster g1 where a is; if the third and second lines are 2 410, both reads (a=2, b=4) are recorded, then all reads in the cluster where b is located are classified into the cluster where a is located, and finally the clustering of all HiFi Reads is completed.
[0081] In the final result acquisition module: In some embodiments, based on the clustering results of high-precision long reads, a certain type of high-precision long reads sequencing data in the clustering results is used as input, and the sequencing data of this type of clustering results are assembled using HiFiasm software, and then the high-precision long reads of other clustering results are assembled using the same method, and finally the results of merging different clustering results are merged to obtain the final assembly result.
[0082] For example, g1: 1, 2, 3, g2: 4, 5, 6; based on the read IDs contained in g1 and g2, g1.fa and g2.fa are extracted from the total reads, representing the sequencing data of different clustering results, respectively. g1.fa and g2.fa are used as input, and HiFiasm software is used to assemble g1.fa and g2.fa separately. Then, g1.fa and g2.fa are merged into gcombine.fa as the final result of the divide-and-conquer assembly method for different chromosomes.
[0083] In a third aspect, the present application provides a computer device comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the genome assembly method based on the Perl language as described above.
[0084] The computer device may be a desktop computer, a notebook computer, a PDA, a cloud server, etc. The computer device may interact with the user via a keyboard, a mouse, a remote control, a touchpad, or a voice control device.
[0085] The memory includes at least one type of readable storage medium, including flash memory, hard disk, multimedia card, card-type memory (e.g., SD or D interface display memory, etc.), random access memory (RAM), static random access memory (SRAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), programmable read-only memory (PROM), magnetic storage, magnetic disk, optical disk, etc. In some embodiments, the memory can be an internal storage unit of the computer device, such as the hard disk or memory of the computer device. In other embodiments, the memory can also be an external storage device of the computer device, such as a plug-in hard disk equipped on the computer device, a smart memory card (SmartMediaCard, SMC), a secure digital (SecureDigital, SD) card, a flash memory card (FlashCard), etc. Of course, the memory can also include both the internal storage unit of the computer device and its external storage device. In this embodiment, the memory is commonly used to store the operating system and various application software installed on the computer device, such as the program code of the genome assembly method based on the Perl language, etc. In addition, the memory can also be used to temporarily store various types of data that have been output or are about to be output.
[0086] The processor can be a central processing unit (CPU), a controller, a microcontroller, a microprocessor or other data processing chips in some embodiments. The processor is generally used for controlling the overall operation of the computer device. In the present embodiment, the processor is used for running the program code or the processing data stored in the memory, for example, running the program code of the genome assembly method based on the perl language.
[0087] In a fourth aspect, the present application provides a computer-readable storage medium storing a computer program, which, when executed by a processor, performs the steps of the genome assembly method based on the Perl language.
[0088] The computer-readable storage medium stores an interface display program, and the interface display program can be executed by at least one processor to enable the at least one processor to perform the steps of the genome assembly method based on the Perl language as described above.
[0089] The above is a detailed introduction to the genome assembly method based on the Perl language provided in the examples. Specific examples are used herein to illustrate the principles and implementation methods of the present application. The description of the above examples is only used to help understand the method of the present application and its core idea. At the same time, for those skilled in the art, according to the ideas of the present application, there may be changes in the specific implementation methods and application scopes. In summary, the content of this specification should not be understood as limiting the present application.
Claims
1. A genome assembly method based on the Perl language, characterized by: The steps include: Obtain a high-precision short Reads set, and split the 1 to SRL sites of each Read in the high-precision short Reads set into SRL-L+1 small fragments of length L, which are recorded as S set, where the elements are S i ; By S i The distribution of depth determines the average sequencing depth AD of the high-precision short read set, and the sequencing depth is calculated by the interval from AD to the inflection point that does not conform to the Poisson distribution data, which conforms to the standard deviation STD of the Poisson distribution data; A normal distribution model was used to fit the theoretical expected value of sequencing depth based on AD and STD. The upper limit Emin and lower limit Emax of the theoretical expected sequencing depth were defined based on the coverage area reaching 95% on both sides of the normal distribution mean. Based on the size of Emin and Emax, Si in the range from the upper limit Emin to the lower limit Emax is filtered out from the S set, which is recorded as the filtered Sf set; Get the first effective Reff i ; Get a high-precision long Reads set and record each Read in the high-precision long Reads set as HR i , based on the first effective Reff i Get the second effective Reffi and determine the second effective Reff i The corresponding high-precision long Reads in the genome, if the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN; According to the level of Reff_set_filterN, keep the same as HR i Corresponding to the high-precision long Reads with the largest number of Sfs, the high-precision long Reads corresponding to Reff_set_filter are recorded as Reff_set_filter_best; According to Reff_set_filter_best, get any two HRs i The original interaction strength two-dimensional matrix ForceMatrix; Use any value v in ForceMatrix ij The lengths of the two corresponding high-precision long reads are corrected to obtain the corrected high-precision long read interaction strength two-dimensional matrix NormalForceMatrix; Split the two-dimensional matrix NormalForceMatrix to obtain high-precision clustering results of long reads; According to the clustering results of high-precision long reads, assembly is performed to obtain the final results.
2. The genome assembly method based on the Perl language according to claim 1, characterized in that: Obtain a Hi-C sequencing Reads set, and split the 1 to HRL sites of each Read in the Hi-C sequencing Reads set into HRL-L+1 subsets to obtain the Sf of each Read. i Set, the number of sets is NSf i , Sf i The number of occurrences in the Sf set is recorded as ENSf i , use ENSf i Divide by NSf i , and obtain the ratio Re. When Re is greater than or equal to 0.8, it is considered that the read in the Hi-C sequencing reads set is suitable for identifying the clustering relationship between different reads in the high-precision long reads set, and is recorded as the first effective Reff i and / or Get a high-precision long Reads set and record each Read in the high-precision long Reads set as HR i , HR i 1 to LRL sites are divided into LRL-L+1 small segments of length L and recorded as set HRSET i and / or According to the different HRs corresponding to the forward and reverse reads in Reff_set_filter_best i Count the cases, record them as Hai, Hbi, and iterate them pair by pair to get any two HRs. i The two-dimensional matrix ForceMatrix of the raw interaction strength between them; and / or The two high-precision long reads are recorded as Hifi1 and Hifi2; the lengths of the two high-precision long reads are recorded as LengthHifi1 and LengthHifi2, and the correction formula is v ij *1e8 / (LengthHifi1+LengthHifi2), according to the correction formula, the corrected high-precision long reads interaction strength two-dimensional matrix NormalForceMatrix is obtained.
3. The genome assembly method based on the Perl language according to claim 2, characterized in that: Obtain a high-precision short read collection of the target species through high-throughput sequencing; and / or Obtaining a Hi-C sequencing read set of the target species through high-throughput sequencing; and / or High-throughput sequencing is used to obtain a high-precision long read collection of the target species.
4. The genome assembly method based on the Perl language according to claim 2, characterized in that: The first valid Reff is recorded by the hash data structure Hash-Reff i The number of elements in the Sf set ENSf i , ENSf i Sort from large to small, and the top 95% of the values are recorded as the second effective Reffi. The second effective Reffi is calculated in turn with the intersection of HRSTi to obtain the number of elements Ninter contained in the intersection. If the number Ninter is less than half of the second Reffi length HRL, it is recorded as 0, and the rest are counted according to the actual number. If the high-precision long reads corresponding to the forward and reverse Hi-C sequencing reads come from the same line, then the second effective Reffi is considered i It does not provide effective information for clustering high-precision long Reads; if not, the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN.
5. The genome assembly method based on the Perl language according to claim 2, characterized in that: In the pair-by-pair iteration process, if HR i If it is greater than the existing Hai, Hbi, the different HRs corresponding to the forward and reverse reads will be i replace.
6. The genome assembly method based on the Perl language according to claim 2, characterized in that: The two-dimensional matrix NormalForceMatrix is split into a list of corrected intensities between any two high-precision long reads. The first and second columns represent two high-precision long reads a and b, and the third column represents the interaction strength v. The clustering process first checks v and obtains a descending order of v. When the v value is the largest, the corresponding two high-precision long reads are first picked out to create a cluster and recorded. Then, v is checked in sequence. If the Reads ID in the row has been recorded, it is included in the clustered group. If it has not been recorded, a new cluster is created and recorded again to finally obtain the clustering result of the high-precision long reads.
7. The genome assembly method based on Perl language according to claim 2, characterized in that: According to the clustering results of high-precision long reads, the sequencing data of a certain type of high-precision long reads in the clustering results were used as input, and the sequencing data of this type of clustering results were assembled using HiFiasm software. Then, the high-precision long reads of other clustering results were assembled using the same method. Finally, the results of merging different clustering results were merged to obtain the final assembly result.
8. A genome assembly system based on the Perl language, characterized in that: include: The high-precision short reads collection module is used to obtain a high-precision short reads collection. The 1 to SRL sites of each read in the high-precision short reads collection are divided into SRL-L+1 small fragments of length L, which are recorded as S sets. The elements of the S i ; AD and STD modules for S i The distribution of depth determines the average sequencing depth AD of the high-precision short read set, and the sequencing depth is calculated by the interval from AD to the inflection point that does not conform to the Poisson distribution data, which conforms to the standard deviation STD of the Poisson distribution data; The upper limit Emin and lower limit Emax position determination module is used to fit the theoretical expected value of sequencing depth based on AD and STD using a normal distribution model, and define the upper limit Emin and lower limit Emax of the theoretical expectation of sequencing depth based on the region where the coverage area on both sides of the normal distribution mean reaches 95%; The module for obtaining the filtered Sf set is used to filter out Si within the range of the upper limit Emin to the lower limit Emax from the S set based on the size of Emin and Emax, and record it as the filtered Sf set; First effective Reff i Get module, used to get the first effective Reff i ; Reff_set_filterN is a module used to obtain a high-precision long Reads set, and each Read in the high-precision long Reads set is recorded as HR i , based on the first effective Reff i Get the second effective Reffi and determine the second effective Reff i The corresponding high-precision long Reads in the genome, if the second effective Reff i For high-precision long Read clustering, it can provide effective information, recorded as Reff_set_filter, and record the second effective Reff i The number of overlapping Sf between different high-precision long reads Reff_set_filterN; Reads is the module for obtaining Reff_set_filter_best, which is used to retain the value of HR according to the level of Reff_set_filterN. i Corresponding to the high-precision long Reads with the largest number of Sfs, the high-precision long Reads corresponding to Reff_set_filter are recorded as Reff_set_filter_best; ForceMatrix acquisition module, used to obtain any two HRs according to Reff_set_filter_best i The original interaction strength two-dimensional matrix ForceMatrix; NormalForceMatrix obtains the module, using any value v in ForceMatrix ij The lengths of the two corresponding high-precision long reads are corrected to obtain the corrected high-precision long read interaction strength two-dimensional matrix NormalForceMatrix; The module for obtaining the clustering results of high-precision long reads is used to split the two-dimensional matrix NormalForceMatrix to obtain the clustering results of high-precision long reads; The final result acquisition module is used to assemble the final results based on the clustering results of high-precision long reads.
9. Computer device, characterized in that The device comprises a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the method according to any one of claims 1 to 7.
10. A computer-readable storage medium, characterized in that A computer program is stored, and when the computer program is executed by a processor, the processor is caused to perform the steps of the method according to any one of claims 1 to 7.