JP2005218421A - Method for producing contig from dna fragmented sequential data by total genom shot gun method and recording medium - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a practical contig by classifying correct annexations, incorrect annexations and annexations which can not be judged as correct or incorrect substantially, from fragmented sequential data obtained by a total genom shot gun method without using already known repeated sequential data in a genom sequential determination containing repeated sequences. <P>SOLUTION: This method for producing the contig is provided by examining whether each of the fragmented sequences is included in another fragmented sequence approximately, and in performing the annexation on a fragmented sequence x not included in any other fragmented sequences and its complemental sequence, determining a fragment sequence or its part capable of extending the fragmented sequence x in 3' direction based on a score calculated from a partial fragmented sequence which is in common with the fragmented sequence x, among another fragmented sequence capable of extending the fragmented sequence in 3' direction or its complemental sequence. <P>COPYRIGHT: (C)2005,JPO&NCIPI

Description

  The present invention relates to the field of genome base sequencing methods, and more particularly to a method and a recording medium for creating a contig from a data set including DNA fragment sequences by the whole genome shotgun method.

  The whole genome shotgun method is widely used for genome sequencing. In the whole shotgun method, the genome is first physically divided into fine fragments, the obtained fragments are randomly cloned, and the nucleotide sequence of each obtained clone is determined. The fragment sequence thus obtained contains a nucleotide sequence at a random position in the genome sequence. Next, a sequence shared by the fragment sequences is detected, the fragment sequences sharing the sequences are merged to create a contig, and the original genomic base sequence or a part thereof is restored.

  However, the genome sequence may contain a large number of repetitive sequences, and inconsistencies may occur when fragment sequences having a common sequence are merged randomly. Therefore, there is a need for a method for merging fragment sequences so that no contradiction occurs even if there are repetitive sequences in the genome sequence.

  In Non-Patent Document 1, Myers et al. Erroneously masked a repetitive sequence in a fragment sequence having homology with a known repetitive sequence using a database of known repetitive sequences and did not align the masked region. A method is proposed to avoid sequence merging. In this case, a boundary between a repetitive sequence and a non-repetitive sequence, that is, a repetitive sequence boundary is detected, and a method for prohibiting the merging of fragment sequences across the boundary is used together, and fragment sequences that do not cause contradiction are merged.

  The group of Batzoglou et al. In Non-Patent Document 2 has devised a method for reducing erroneous merging due to repetitive sequences by preferentially merging mate pairs that overlap in pairs with detection of repetitive sequence boundaries.

  In Non-Patent Document 3, the group of Wang et al. Is derived from a repetitive sequence with respect to a k-mer having an appearance frequency of a certain level or more by determining the appearance frequency of a k-mer (k character base sequence) appearing in a fragment sequence. In view of this, a method has been proposed in which a k-mer that is derived from a repetitive sequence is masked in a fragment sequence, and an alignment score is not given to the masked region, thereby reducing erroneous merging.

Eugene W. Myers, Granger G., et al. Sutton et al, A Whole Genome Assembly of Drosophila, Science (2000) (287) 2196-2204 Serafim Batzoglou et al, ARACHNE: A Whole-Genome Shotgun Assembler, Genome Research, January 2002 (1): pp. 177-189 Jun Wang et al, A Sequence Assembler That Masks Exact Repeats Isolated from the Shotgun Data, Genome Research, March 2002 (page 8): 82.

  However, the method of masking repetitive sequences cannot be applied unless the repetitive sequences are known. For an organism for which a new genomic sequence is to be determined, the repetitive sequences present in that organism are unknown, and there is often no database for known repetitive sequences. In addition, even in a species for which a database for known repetitive sequences is maintained, the repetitive sequence database is created based on the result of determining a part of the genome sequence. Normally, only repetitive sequences are included in the database, and repetitive sequences not listed in the database may cause incorrect fragment merging.

  The method of masking repetitive sequences using the frequency of k-mer appearing in the fragment sequence does not require a repetitive sequence database. However, as described in Non-Patent Document 3, it is possible to mask a repetitive sequence having a very high number of occurrences in the genome with a high probability, but erroneously determining a repetitive sequence having a relatively low number of appearances. Often, it would cause incorrect fragment merging.

  In addition, the method for detecting repetitive sequence boundaries has been difficult to apply to genome sequences in which repetitive sequences are frequently scattered. This is because repeated sequence boundaries appear frequently, so that fragment sequences are hardly merged, resulting in a very large number of contigs and impractical. In addition, even the repetitive sequences appearing in the genome at a low frequency are sometimes suppressed in merging which is essentially possible, which increases the number of contigs.

  In the case of Non-Patent Document 1 Myers, the repetitive sequence that is frequently scattered is masked by the method of masking the repetitive sequence, and the repetitive sequence boundary is detected after the repetitive sequence remaining without being masked is suppressed to a low frequency. This is handled by using a method that requires a repetitive sequence database.

  The group of Batzoglou et al. In Non-Patent Document 2 reduces the influence of high-frequency repetitive sequences by preferentially merging mate pairs that overlap in pairs, but as described in Non-Patent Document 2 The effect is considered to be limited, such as the possibility of erroneous merging when one of the pair of mate pairs does not contain a non-repetitive sequence.

  Because higher organisms often have genomic sequences that contain a high number of scattered repeats, determining the genome sequence for these organisms does not require a database of known repeats, There has been a demand for creating a contig without erroneous merging even if repeated sequences scattered in the frequency are present in the genome sequence.

  Based on a data set containing base sequence information of random DNA fragments acquired by the whole genome shotgun method, it is checked whether each fragment sequence is approximately included in other fragment sequences, and included in any fragment sequence A fragment sequence x and its complementary sequence that are not merged, and a portion shared with the fragment sequence x from another fragment sequence or its complementary sequence that can extend the fragment sequence in the 3 ′ direction when merged The problem is solved by a method for determining a fragment sequence or a portion thereof that extends the fragment sequence x in the 3 ′ direction based on the alignment score calculated from the fragment sequence.

  More preferably, a score proportional to the length of the homologous portion is used instead of the alignment score.

  More preferably, with respect to each fragment sequence x, the processing speed is improved by a method of calculating the score in order from the fragment sequence having a long potential shared sequence length among the fragment sequences extending the fragment sequence x in the 3 ′ direction. be able to.

  More preferably, data having a hash table in which both continuous or non-contiguous partial base sequences appearing in a fragment sequence and the position of the partial base sequence from the beginning of the fragment sequence are simultaneously associated with a set of fragment sequences is recorded. By using a computer-readable recording medium, it is possible to reduce the time required to obtain a fragment sequence having a necessary shared sequence length.

  The present invention relates to a) merging errors due to repetitive sequences regardless of the presence or absence of a database of known repetitive sequences and the magnitude of the number of occurrences of repetitive sequences in the genome and the distribution of repetitive sequences in the genome. B) errors C) It is possible to classify the three types of merge that cannot be specified essentially from the fragment sequence information, thereby generating a contig with few errors.

  Hereinafter, an embodiment of the present invention will be described with reference to the accompanying drawings.

[Definition]
“Read” means a short sequence of DNA sequence information derived from a random position from the genome. The lead sequence is a sequence having a direction from the 5 ′ end to the 3 ′ end, and the 5 ′ side is the head.

  “Complementary leads” or “two complementary leads” are in a “complementary” relationship means that they are in a DNA double-stranded relationship.

  A unique “lead number” is assigned to all the leads used in this embodiment. It is desirable that leads in a complementary relationship can be uniquely identified by some method. In this embodiment, as one method, when the difference in lead numbers is equal to the number of input leads, the two leads are complementary to each other. The lead number was determined as follows.

  A continuous partial base sequence or a non-continuous partial base sequence in a base sequence is called a “seed”.

  Two leads “duplicate” means that the two lead sequences have a partial sequence that aligns, and the aligned partial sequences terminate at the end of either lead. The partial sequences to be aligned do not need to be completely coincident, and it is determined whether or not the alignment is made according to the parameter of the degree of difference given for each process.

  If two leads are x and y and the two leads x and y overlap, the following five states are considered.

  As shown in FIG. 2, the 5 ′ end of y is on the inside (3 ′ side) from the 5 ′ end of x, and the 3 ′ end of y is on the outside (3 ′ side) from the 3 ′ end of x. In this state, it is called “the overlapping destination of x (in the 3 ′ direction) is y”.

As shown in FIG. 3, the 5 ′ end of y is outside (5 ′ side) from the 5 ′ end of x, and the 3 ′ end of y is inside (5 ′ side) from the 3 ′ end of x. Respective complementary read ^x c of the state x and y, when viewed from the ^{y c,} the 5 'terminus of ^{y c} from the 5' end of the ^{x c} is 'located in (side, 3 ^{x c} inside 3)' from end Since the 3 ′ end of y ^c is on the outside (3 ′ side), it is referred to as “the overlapping destination of x ^c is y ^c ”.

  As shown in FIG. 4, the 5 ′ end of y is inside (3 ′ side) from the 5 ′ end of x, and the 3 ′ end of y is inside (5 ′ side) from the 3 ′ end of x. At this time, a case where only one of the 5 'end 3' ends of the two leads is also included in this state. In this state, “x contains y” or “y is contained in z”, x is called “parent lead” of y, and y is called “child lead” of x. That is, the parent lead includes a child lead.

  As shown in FIG. 5, the 5 'end of y is outside (5' side) from the 5 'end of x, and the 3' end of y is outside (3 'side) from the 3' end of x. At this time, a case where only one of the 5 'end 3' ends of the two leads is also included in this state. In this state, “y contains x” or “x is contained in y”, y is called “parent lead” of x, and x is called “child lead” of y.

  As shown in FIG. 6, the 5 ′ end of x and the 5 ′ end of y are matched, and the 3 ′ end of x and the 3 ′ end of y are also matched. Such a state is also a kind of “include” described above, and the lead having the smaller lead number is the parent lead and the lead having the larger lead number is the child lead.

As for the above inclusion relationship, the same relationship generally holds for each complementary lead. When “x includes y”, “x ^c also includes y ^c ”. However, depending on how the lead numbers are assigned between the complementary leads, this may not necessarily be the case when the leads as shown in FIG. 6 match. However, the relationship between the lead numbers assigned in the present embodiment is established. ing.

[Build hash table]
A lead number is assigned to each lead acquired by the whole genome shotgun method. For each lead, a set of lead numbers of the lead sequence that holds the partial base sequence at the above position, using the seed sequence in the lead sequence and the position from the beginning of the lead sequence of the seed sequence as an index at the same time, Build a hash table that can be obtained. A simplified diagram for illustration is shown in FIG. At this time, the position is counted from the head (5 ′ side) of the lead. The length of the seed sequence is preferably 6 to 14 bases. More preferably, the hash table is implemented using an array structure.

[Exclude frequently repeated sequences]
For sequences that are expected to appear universally at a high frequency in the genome sequence, a read number is not added to the hash table. As one example, for example, a seed sequence such as “AAAAAAAAAA” or “ATATATATAT” such as a simple repeat sequence consisting of 1 base, 2 bases and 3 bases may be excluded from the hash table. However, other implementation methods such as actually determining the frequency at which each partial base sequence appears in the fragment sequence and excluding the lead number from the hash table for the partial base sequence appearing at a certain frequency or higher are also conceivable.

  First, select one lead to be inspected and set it as a. This test sequence is repeated across all reads. The variable w is initialized to 1, and a candidate for a read that is shifted by w bases in the 3 'direction with respect to a and becomes the duplication destination of a is searched from the hash table. Specifically, for a seed sequence existing after the (w + 1) th base in a, a set of read lead numbers such that the seed sequence appears at a position shifted in the 5 ′ direction by w from the position appearing in a. Ask.

[Refinement of alignment with seed as core]
Let b be the lead indicated by each obtained lead number. As shown in FIG. 8, the alignment is refined in both 5 ′ and 3 ′ directions with the seed group shared between a and b as the nucleus. This step is preferably performed by a dynamic programming method such as Smith-Waterman method or its extension, but may be performed by a sparse dynamic programming method or the like.

[Refine processing by multiple seeds]
One simple way to perform the above steps is to use a single seed as the core of the alignment. In the case of such a method, if a seed sequence as a nucleus is selected in a repetitive region in a genome sequence, a processing for a set of reads that does not actually have a significant alignment will be performed. . In order to avoid this and increase the efficiency of processing, it is also conceivable to combine a plurality of seeds instead of a single seed. One way to combine these seeds is that the difference in position coordinates where seeds appear in each lead can be thought of as approximating the deviation of the overlap between the leads. As indicated by 9, it is also conceivable to combine the seeds according to the difference in the appearance position of each lead.

  Depending on how b is selected in the above-described steps, in many cases, the relationship between a and b is the relationship that “a overlaps with b” or the relationship “a includes b”. Sometimes.

[Destination evaluation]
Assume that a lead x (in the 3 ′ direction) overlaps a lead y. There can be a plurality of such leads y for one x. There may be none. For such a plurality of duplication destinations, a scale for evaluating each “goodness as duplication destination” is provided and calculated. As such a measure, the shared sequence length of the overlap of leads x and y, the alignment score of x and y, and the like can be considered.

[Select the best duplication destination]
When the duplication destination b for a certain lead a is obtained by the above steps, if the duplication destination of the lead a is only b, the lead b is set as the best duplication destination at the present time of the lead a. When the duplication destination b for a is obtained, if the best duplication destination c so far already exists, the duplication between a and b is compared with the duplication between a and c according to the above-mentioned scale. The better duplication is set as the best duplication destination of a at this time.

Exclude included leads
If refinement of the alignment results in a situation where “lead a includes lead b” or “lead b includes lead a”, the “included” lead is “included” in the following steps. Since the lead duplication destination is not calculated and is not directly adopted as the duplication destination, information representing this is added to the “included” lead.

[Potential shared sequence length]
In addition, as shown in FIG. 10, in the above method, a read number that is approximately predicted to have an alignment that is shifted by about w bases from the beginning of the lead a is obtained from the hash table. Is expected to have a shared sequence approximately “a sequence length−w” in length. This is called the “potential shared sequence length” of the two reads that align.

[Estimation of good overlap obtained]
By starting w from 1, processing is performed in order from a read group having a long potential shared sequence length for a. Therefore, it is determined whether the duplication destination better than the current best duplication destination can be obtained by continuing this process. For example, if the alignment shared sequence length is used as a measure of overlap, the potential shared sequence length can be considered an approximation of the exact shared sequence length of this alignment, and as the process continues Since the potential shared sequence length of the lead to be processed with the lead a is monotonously shortened, even if processing is continued thereafter, a better duplication destination than the best duplication destination at the present time cannot be obtained. It can be judged. Specifically, potential shared sequence length × (1 + error tolerance) × alignment match score or 1 (when the shared sequence length is adopted in the scale) is below the score for the current best duplication destination. It is determined that there is no better duplication destination. For the error tolerance, an estimate of the upper limit of the proportion of polymorphism expected in the target genome is given in advance.

[End iteration]
If it is determined that even if processing is continued further than this, it is not possible to obtain a better duplication destination than the current best duplication destination, the processing loop for w is terminated. Thereby, unnecessary calculation processing can be omitted.

  When the duplication destination has not been obtained yet or when it is determined that there is a possibility that a better duplication destination may be obtained by continuing the processing, w is increased and the above steps are repeated. At this time, the increment of w is preferably the sequence length of the seed used in the construction of the hash table, but it is conceivable that the increment is changed according to the best overlap alignment obtained at present.

  If it is determined that the lead a itself is an “included lead”, there is no need to obtain the duplication destination, so the step for a is immediately terminated, and the next search target sequence is set as a new a. Iterate.

[True parent lead, true child lead]
A lead that has not been determined to be an “included lead” by the above steps is a lead that has never been determined to be a “child lead”. Such a lead is called a “true parent lead”. A lead that is determined to be a “child lead” at least once is called a “true child lead”.

[Calculation of inclusion relationship]
Next, a process for calculating a true parent lead including the lead is performed for those determined to be true child leads in the above steps.

  This step performs almost the same process as the above-described step of searching for duplication destinations, but the inspection target array is already a true child read, and the alignment target indicated by the set of lead numbers obtained by searching the hash table The only lead is the true parent lead.

  In this step, the data structure indicating the relationship between the true parent lead and the true child lead is calculated by the same process as in the previous step, but in this step, one “true parent lead that includes me” was found. Do not end the iteration at that point and continue until you find all the “true parent leads that contain you”.

  By this step, if a certain true child lead is c, c has information of one or more “true parent lead of c”, and if a certain true parent lead is p, p is “true child lead”. Do not have "or [one or more true child leads] information.

  As described above, the duplicate destination of the true parent lead is determined for each true parent lead when it exists.

[Handling duplicates to true child leads]
As shown in FIG. 11, when there is no duplication destination of the true parent lead for each lead that is a true parent lead, but the true child lead exists as the duplication destination, alignment is performed as shown in FIG. The parameters of the elaboration process are relaxed, and the step of calculating the overlap with the true parent lead that is aligned from the middle to the test target sequence is performed, and the true overlap with the best overlap among such alignment is performed. The portion of the parent lead that is shared with the lead to be inspected is used as an overlapping destination array.

  As a result of the above processing, one or zero 3′-direction duplication destination sequences are determined for all true parent reads and their complementary sequences as shown in FIG. This can be viewed as a directed graph with each array as a point and a relationship indicating an overlap destination in the 3 'direction as an edge.

[Graph integration for complementary sequences]
As shown in FIG. 1, by integrating the information of complementary sequences for each sequence, a graph showing the overlapping relationship in both the 3 ′ direction and the 5 ′ direction of each base sequence is derived. In this integration, when the overlapping destination in the 3 ′ direction of a certain sequence j is the sequence k and the overlapping destination in the 3 ′ direction of the complementary sequence of k is the complementary sequence of j, the edge connecting j and k Are merged into one.

[Calculation of number of references]
For the graph constructed by the above processing, it is calculated how many sides each point corresponding to the base sequence holds for each of the 3 ′ direction and the 5 ′ direction. A point having two or more sides in at least one of the 3 ′ direction and the 5 ′ direction is called a “branch point”.

[3 'direction merge inspection]
As shown in FIG. 13, in the graph after integration, an arbitrary point that has not been reached and is not a branch point in the following steps is used as a starting point, the side in the 3 ′ direction is traced, and a portion that is not a branch point is Consecutive pairs.

[5 'direction merge inspection]
Further, as shown in FIG. 14, the same processing as described above is performed in the 5 ′ direction. At this time, with respect to the branch point, whether to merge with this processing set is determined by the number of the respective sides in the 5 ′ direction and 3 ′ direction of the branch point. In the case of processing from a direction having two or more sides when viewed from the branch point, the branch point is not included in the set.

[Restore True Child Lead]
Through the above steps, the true parent lead or its partial arrangement can be divided into sets for each part having no branch on the graph. For short repetitive sequences, a graph as shown in FIG. 16 is generated, and the upper and lower two types of read groups are never mixed. For a long repetitive sequence, a graph as shown in FIG. 17 is generated, and two types of upper and lower leads are connected through a branch point. When a branch point occurs, merging cannot be performed at the position of the branch point. For example, in the example of FIG. 17, even if the sequence A and the sequence C are reversed, the fragment sequences obtained by the shotgun method may be the same, resulting in ambiguity. The set of arrays obtained by the above steps constitutes a column of graph edges and overlapping relationships. For a lead that has been determined to be a true child lead in the previous processing for this array column, the true parent lead is single, or all of the true parent leads are classified into the same set and are adjacent. Only in the case of an array, it is added to this array set as shown in FIG.

  The set of sequences separated by the above steps can be aligned from the 5 ′ direction to the 3 ′ direction by using the information representing the position they have, and the contig sequence can be obtained by multiple alignment. Generate.

It is a figure showing the graph structure which shows the branch of a lead. It is a figure showing the state which overlaps with lead y from lead x. It is a figure showing the state which overlaps with lead x from lead y. It is a figure showing the state where lead x includes lead y. It is a figure showing the state where lead y includes lead x. It is a figure showing the state where lead x is equal to lead y. It is a figure showing the hash table which makes an index the seed arrangement | sequence and the position in an arrangement | sequence. It is a figure showing refinement of alignment by the nucleus of lead a and lead b. It is a figure showing the classification | category of a seed group. It is a figure showing potential heavy rain arrangement length. It is a figure showing the detection of the duplication state to the true child lead and the duplication to the true parent lead. It is a figure showing the graph of the duplication relation of a lead. It is a figure showing the process to 3 'direction from a starting point. It is a figure showing the process to 5 'direction from a starting point. It is a figure showing restoration of a true child lead. An example where the correct duplication destination over a short repeat sequence is determined. An example with long repeat sequences and inherent ambiguity.

Claims (5)

  (A) Check whether each fragment sequence is approximately included in other fragment sequences against a data set including base sequence information of random DNA fragments acquired by the whole genome shotgun method (b) which fragment For a fragment sequence x not included in the sequence and its complementary sequence, another fragment sequence that can extend the fragment sequence in the 3 ′ direction when merged or a complementary sequence thereof, between the fragment sequence x A method for determining a fragment sequence or a portion thereof that extends a fragment sequence x in the 3 ′ direction based on an alignment score calculated from a shared partial fragment sequence.   The method according to claim 1, wherein a score proportional to the length of the homologous portion is used instead of the alignment score of (b).   The method according to claim 1 or 2, wherein for each fragment sequence x, the score is calculated in order from a fragment sequence having a longer potential shared sequence length among fragment sequences extending the fragment sequence x in the 3 'direction. How to do the processing.   Computer-readable data recorded with a hash table that associates both a continuous or non-contiguous partial base sequence appearing in a fragment sequence with the position of the partial base sequence from the beginning of the fragment sequence as an index simultaneously with a set of fragment sequences Recording medium.   A computer-readable recording medium recording a program for causing a computer to execute the method according to claim 1.

Images