In bioinformatics, scoring matrix is also called as substitution matrix. As for protein sequences, the matrix is often based on observed substitution rates, derived from the substitution frequencies seen in multiple alignments of sequences. Every possible identity and substitution is assigned a score based on the observed frequencies in alignments of related proteins. Similarly, every possible identity and substitution in alignments of related DNA sequences is also assigned a score. However, the scoring matrix in alignment of DNA sequences is relatively simple and intuitional. Table 1 is a scoring matrix. An identity in scoring matrix is assigned a positive score R (R>0). A substitution is assigned a positive score S (S>0). But the score S must be subtracted from the total score of an alignment. For example, an alignment of two short DNA sequences ATGGTGCAACTGACT and ATGGTGCACTTGACT is the following:The score of alignment should be 13R−2S.Table 1

A usual scoring matrix.

Another important problem for alignment is the treatment of gaps, i.e., spaces inserted to optimize the alignment score. A ‘gap open’ penalty is one that is the cost for the first space of each gap spaces. A ‘gap extension’ penalty is one that is the cost for one of each gap spaces except for the first space. Typically, the cost of extending a gap is set to be 5–10 times lower than the cost for opening a gap (http://www.ebi.ac.uk/Tools/emboss/align/index.html). There is one way to compute a penalty for a gap of n positions: gap opening penalty + (n−1)* gap extension penalty. Now, let parameter a be gap opening penalty and parameter b be gap extension penalty. And let parameter k be n−1. Then we have the scoring equation: R0*R−S0*S−T0*(a+bk), where R is the score of each match, and S is the score of each mismatch and a+bk is the score of each gap. The parameters R ₀, S ₀ and T ₀ denote the total amount of matches, the total amount of mismatches and the total amount of gaps, respectively.

There are two primary DNA sequences: G ₁ (g1g2⋯gM) and G ₂ (g1g2⋯gN), where M and N denote the length of G ₁ and G ₂, respectively. Fig. 1 shows the building steps of all possible simple alignments without spaces within G ₁ and G ₂. The gap formed by spaces lies in the hanging ends of the overlap. Step 1 of Fig. 1 gives the initial alignment that the last base of G ₁ overlaps the first base of G ₂. Then every time G ₁ moves one base position along the direction of G ₂. Step M+N−1 of Fig. 1 gives the last alignment that the first base of G ₁ overlaps the last base of G ₂. Every alignment has a score according to the scoring equation R0*R−S0*S−T0*(a+bk). Then we can obtain a serial of dots (x,y), where x denotes index of steps and y denotes the score value corresponding to x. When one connects adjacent dots with lines, then one obtains a zigzag like curve of definite geometrical shape. In Fig. 2 we illustrate the graphical score representation of simple alignments of sequences G ₁ and G ₂ (G ₁, GGCCTCTGCCTAATCACACAGATCTAACAGGATTATTTC; G ₂, GGCCTCT GCCTTATTACACAAATCTTAACAGGACTATTTC). The scoring equation is R0*R−S0*S−T0*(a+bk), where identity score R is 9, substitution score S is 1, gap opening penalty a is 15, and gap extension penalty b is 1.Fig. 1

The building steps of all possible alignments without spaces within G₁ and G₂.

The above simple alignment algorithm shown in Fig. 1 is a time-consuming algorithm. Its time complexity is O(N2) if the two aligned sequences have equal length N. Here, we provide an improved alignment process to obtain a lower time complexity. A close look at Fig. 1 shows that the beginning steps and the ending steps of the sliding process are not necessary. These unnecessary steps consume some time. Now, we need to find out the unnecessary steps to save time.

Fig. 3 shows the improved process. Step (1) of Fig. 3 gives the initial alignment that the first base of G ₁ overlaps the first base of G ₂. Then the sliding process is done along the left and the right, respectively. Firstly, we consider the sliding process is done along the right, as shown in step (2). Every alignment has a score S according to the scoring equation R0*R−S0*S−T0*(a+bk). Let score S_i be the score of the ith step. Let Smax be the top score within all the scores (i.e., Smax=max{S1,S2,…,Si−1}, i>1). As for every step, we still define another score Si′. The score Si′ is obtained by another scoring equation R0*R+S0*R−T0*(a+bk). In the new scoring equation we assume all overlapped bases are matched each other. The sliding process is stopped when Smax≥Si′. Obviously, the remaining steps along the right are unnecessary because all scores of the remaining steps have no chance to obtain a higher score than Smax.Fig. 3

Improved simple alignment process with less consuming time (The length of G₁ is M. The length of G₂ is N. And let M≤N).

By the simple alignment algorithm a good simple alignment R of G ₁ and G ₂ is determined when the score of the simple alignment reaches its peak. Let C be the largest common substrings of R, where C=ϕ∪{C1,C2,…,Cm}. If C=ϕ, there is no largest common substring. When C={C1,C2,…,Cm}, there are m largest common substrings, where |C1|=|C2|=⋯=|Cm|=k and k denotes the number of matches within a largest common substring. Here, we can see that the largest common substring C_i (i=1,2,…,m) of R devotes its maximal score to the initial alignment reaching peak score. The alignment shows the approximate overall alignment feature of G ₁ and G ₂. However, some C_i of the largest common substrings of R may be a part of a larger common substring than C_i because of the limit of the initial alignment reaching peak score. In order to protect a larger common substring than C_i from being split by C_i and find out the substring, an extension algorithm for the largest common substring is suggested as the following.(1)

Let K be the length of the largest common substring of G ₁ and G ₂. Its value is generated when the simple alignment algorithm applies to the sequences G ₁ and G ₂.

Now, we give a practical example to illustrate the above steps for a better understanding. There are two random sequences: G ₁ (GCCTAGTTCCCCCA) and G ₂ (GCCTCGCATCCCCCA). By the simple alignment algorithm a good simple alignment R of G ₁ and G ₂ is determined, where the alignment R is GCCTAGTTCCCCCAGCCTCGCATCCCCCA.The largest common substring C of R is “GCCT” (C ₁) and “CCCC” (C ₂). However, the largest common substring of G ₁ and G ₂ is “TCCCCCA”. Its length K is 7. Obviously, there exists at least a larger common substring than “GCCT” or “CCCC”. As for “GCCT”, the two new sub-sequences are “GCCTAGTTC” (S11) and “GCCTCGCAT” (S12), respectively. Because there is no increment of sore, the original common substring “GCCT” (C ₁) is unchanged. As for “CCCC”, the two new sub-sequences are “AGTTCCCCCA” (S21) and “CGCATCCCCCA” (S22), respectively. Because there is an increment of sore, the original common substring “CCCC” (C ₂) is replaced with the new largest common substring “TCCCCCA”.

In sum, the extension algorithm for largest common substring protects a larger common substring than C_i from being split by C_i.

Now, we let U denote the substrings spaced by C, where U=ϕ∪{U1,U2,…,Un}. Then the strings G ₁ and G ₂ are aligned into two types of substrings: (i) The largest common substrings C with continuous matches: G1[i]=G2[i] (of course, a single match is also permitted) and (ii) The substrings U with continuous mismatches: G1[i]≠G2[i] and both without spaces. Obviously, a good simple alignment R of G ₁ and G ₂ is alternately organized by C and U (e.g. R=C11U21C31U41, where |C11|=|C31|. Here, the superscript of C or U denotes the level of sub-alignment. The superscript “1” denotes the first level. And the subscript denotes the index of substrings in the current sub-alignment. The alignment R is the result of the first level sub-alignment).

In order to expressly explain the above parameters, we give a practical example. There are two random sequences: G ₁ (GCCTAGTTCCCCCA) and G ₂ (GCCTCGCATCCCCCA). The G1[i] and G2[i] denote the base of G ₁ and G ₂, respectively. They become a match when G1[i]=G2[i]. By the simple alignment algorithm a good simple alignment R of G ₁ and G ₂ is determined, where the alignment R is GCCTAGTTCCCCCAGCCTCGCATCCCCCA.And by the extension algorithm the largest common substring C of R is “GCCT” and “TCCCCCA”. Then G ₁ and G ₂ are organized by C into three substrings: C11 is “GCCT”, U21 is “AGT” and “CGCA”, C31 is “TCCCCCA”.

Given sequences G ₁ and G ₂, a graphical simple alignment (GSA) tree is built up by the aforesaid simple alignment algorithm and extension algorithm for largest common substring. In the following, we will show how to use the algorithms to construct a GSA tree of G ₁ and G ₂.

Let R be a good simple alignment of G ₁ and G ₂. Let Ci1 and Uj1 be the largest common substring of R and be the substring spaced by Ci1, respectively. We can see that the largest common substring Ci1 of R devotes its maximal score to the global alignment. However, the Uj1 of R may provide a increment of the score to the global alignment if given appropriate gaps within Uj1, though the scores of these Uj1 are low in the first level sub-alignment.

In order to explore the appropriate gaps in Uj1, we give the following several operation steps.(1)

Compute the scores of all simple alignments of Uj1 by the simple alignment algorithm. A good simple alignment Rj1 of Uj1 is generated when its score reaches its peak.

We repeatedly do the above steps until all U in the last level sub-alignment cannot be further decomposed by the simple alignment algorithm and the extension algorithm. Then we can obtain a graphical simple alignment tree for strings G ₁ and G ₂, consisting of a series of substrings. The Fig. 4 illustrates an example of GSA tree.Fig. 4

An example of graphical simple alignment tree (Strings G₁ and G₂ includes 6 substrings in the first level sub-alignment: U11C21U31C41U51C61. The substring U11 includes 3 substrings in the second level sub-alignment: C12U22C32. The substring U31 includes 2 substrings in the second level sub-alignment: U42C52. The substring U51 includes 3 substrings in the second level sub-alignment: U62C72U82. The substring U22 includes 2 substrings in the third level sub-alignment: U13C23. The substring U82 includes 2 substrings in the third level sub-alignment: C33U43).

We can obtain a global alignment of strings G ₁ and G ₂ when their GSA tree is generated. In the following, we will show how to use GSA tree to generate a global alignment. Observing the GSA tree (e.g. Fig. 4), we can see that there are two types of nodes: inner nodes and leaf nodes. The inner nodes consist of substrings U that can be aligned to more substrings. The leaf nodes include substrings C and U, where U cannot be further aligned by the GSA tree method. The global alignment of strings G ₁ and G ₂ should be composed of all leaf nodes. In order to obtain the global alignment, the GSA tree is traversed by post-order traversal of tree. Then all inner nodes are deleted from the result of post-order traversal. We will achieve the global alignment of strings G ₁ and G ₂. For example, the result of post-order traversal of Fig. 4 is C12U13C23U22C32U11C21U42C52U31C41U62C72C33U43U82U51C61.The global alignment is C12U13C23C32C21U42C52C41U62C72C33U43C61 by removing all of inner nodes.

In this section, we give three examples to see the validity of GSA tree for the global alignment of 2 DNA sequences. As for every example, we compare the results by GSA tree with the results by the important heuristic approach “EMBOSS Pairwise Alignment Algorithms—needle” (http://www.ebi.ac.uk/Tools/emboss/align/index.html).

We firstly apply the proposed method to the discussed 2 sequences: G ₁ (GGCCTCTGCCTAATCACACAGATCTAACAGGATTATTTC) and G ₂ (GGCC TCTGCCTTATTACACAAATCTTAACAGGACTATTTC). The scoring equation is R0*R−S0*S−T0*(a+bk), where identity score R is 9, substitution score S is 1, gap opening penalty a is 15, and gap extension penalty b is 1. Of course, one may also choose other values as the parameters of the scoring equation according to practical requirement. Here, we choose the same values in order to keep the context consistency. In Fig. 2, we have described the graphical score representation of simple alignments of sequences G ₁ and G ₂. One can easily find out that the score of the simple alignment reaches its peak when index of step in Fig. 1 is 39. Then we can obtain the first level sub-alignment results according to the simple alignment algorithm and the extension algorithm. By similar rules, we can get all possible sub-alignment results. Fig. 5 shows the GSA tree to be used to align the sequences G ₁ and G ₂. Table 2 lists all substrings of Fig. 5 and the max score of each pair of substrings. Then we obtain the global alignment of the sequences: C11U14C24U34C23U44C54U64C22U43C53. According to the results of Table 2, we illustrate the global alignment results as the following:The notation “‐” denotes the gap within sequence. Then we apply “needle” program to G ₁ and G ₂. The values about gap opening penalty a and gap extension penalty b are the same as the values proposed the GSA tree algorithm. The default values about identity score R and substitution score S are chosen because of no other choice. The alignment results are shown as the following:Obviously, we get consistent results by two different methods. Moreover, similar results can be found out in Table 1 of Randić et al. (2006).Fig. 5

The graphical simple alignment tree to be used to align the sequences G₁ and G₂ (G₁, GGCCTCTGCCTAATCACACAGATCTAACAGGATTATTTC; G₂, GGCCTCTGCCTTATTACACAAATCTTAACAGGACTATTTC).

The GSA tree algorithm is a gradual algorithm that step by step explores suitable gaps to improve the score of the alignment between 2 DNA sequences. The scheme firstly uses the simple alignment algorithm and the extension algorithm to construct a GSA tree. The substrings denoted by leaf node of GSA tree are looked on as a part of the final global alignment. Obviously, GSA tree method is a heuristic algorithm. There is no analytical formula to determine the ‘best’ gap values. Though GSA tree algorithm is a heuristic method, and could not always get the optimal alignment, the validity and practicality of the results can be ensured.

Firstly, we discuss the influence about the initial simple alignment of GSA tree algorithm. In Section 3.2, we describe in detail the construction steps of GSA tree by the initial simple alignment. As for the overall result, the initial simple alignment is very important. It determines the extent to which the consolidated results close to optimal results. Next, we discuss the initial simple alignment and its optimized features.

By the simple alignment algorithm a good simple alignment R of sequences G ₁ and G ₂ is determined when the score of the simple alignment reaches its peak. Then the largest common substrings C of R is determined, where C=ϕ∪{C1,C2,…,Cm}. When C={C1,C2,…,Cm}, there are m the largest common substrings. Because of the limit of the initial alignment reaching peak score, some C_i of the largest common substrings of R may be a part of a larger common substring than C_i. An extension algorithm for the largest common substring is suggested to protect a larger common substring than C_i from being split by C_i if the larger common substring exists. Then we replace the original C_i with the new larger common substring also called as C_i. Once the largest common substrings C_i is determined, it is looked on as a part of the overall result. Then in an optimal alignment of G ₁ and G ₂ (This hypothetical alignment may be obtained by an absolutely optimized algorithm), the two sequences of C_i by GSA tree might be entirely overlapping in the optimal alignment, or partially overlapping, or one sequence might be entirely isolated from the other. Obviously, for the latter two cases, a great deal of gaps and mismatches may be introduced. So the substring C_i is most likely to appear in the optimal alignment.

The above explanation gives some approximate analysis instead of strict mathematical reasoning, but the conclusion is reasonable from the biological point of view. As the sequences under comparison are protein coding, gaps with lengths other than multiples of three are highly unlikely, whereas the GSA tree algorithm can avoid many single or two-base gaps by the approximate approach. Unlike “EMBOSS Pairwise Alignment Algorithms—needle (global)” or “—water (local)”, the proposed GSA tree algorithm is to find a balance between the global and local. The algorithm is for aligning two sequences over their entire length. The match areas with local superiority are produced in the global context. In bioinformatics, it may be reasonable to assume that in the global context the local areas with closely related sequences should be reflected because of the stability of these areas. As for two possibly related sequences, giving priority to these local areas may be a better choice. The accurate optimization result could not show these local features because of the global optimization. The proposed GSA tree algorithm is for aligning two sequences over their entire length. This works best with closely related sequences. If one uses GSA tree to align very distantly related sequences, it will produce a result but much of the alignment may have little or no biological significance. The three examples of 4.1 give practical proofs for the validity and practicality of the GSA tree algorithm. The Example 1 gives two very similar but very short sequences. The GSA tree achieves the same results as needle. The Example 3 gives two similar and long sequences. There is only one different area:

This example shows that giving priority to these local areas may be a better choice. The Example 2 gives two random sequences. There are two different areas in their results. A close look at the alignments discovers that both of the methods can find out those very similar local areas. The main differences between them lie in the very distantly related areas.

Finally, the following analysis shows why we consider the substring C_i in the only initial alignment reaching peak score. Because of the limit of the initial alignment reaching peak score, some C_i of the largest common substrings of R may be a part of a larger common substring than C_i. Then an extension algorithm for the largest common substring is suggested to protect a larger common substring than C_i from being split by C_i if the larger common substring exists. In the extension algorithm we introduce three parameters: K, LL and LR, where K be the length of the largest common substring of G1 and G2. When k=K(k=|Ci|), none of the largest common substrings C_i of R can be extended into larger common substring. The largest common substrings of R are C1,C2,…,Cm, respectively. There is no larger common substring than C_i split by C_i. Otherwise, when k<K, there exists at least a larger common substring than C_i. The larger common substring and C_i might be partially overlapping, or one might be entirely isolated from the other. When they are entirely isolated from each other, the larger common substring can be completely preserved and found out in the next sub-alignment. As for being partially overlapping, we face a choice: the larger common substring or C_i. In order to give a better choice, we use parameters K, LL and LR to construct two new sub-sequences Si1 and Si2. We will replace the original C_i with the new substring also called as C_i if there is an increment of sore when the new larger common substring comes into being. So we can protect a larger common substring than Ci1 from being split by Ci1 and eliminate the limitations caused by the only initial alignment reaching peak score.