Shen S., Tuszynski J.A. Theory and Mathematical Methods for Bioformatics
Подождите немного. Документ загружается.
108 3 Modulus Structure Theory
Exercise 15. For the pair of sequences
A = 110102021010322011101302110103210311023232301 ,
B = 1101033320210103220013021100010321022311022301 ,
answer the following questions:
1. If B is considered to be the sequence mutated from A, compute the three
different modes of type-II, type-III, and type-IV mutations, and write
down the corresponding shifting functions.
2. Find the optimal alignment (A
,B
)of(A, B) using the dynamic program-
ming-based algorithm.
3. Give the alignment modulus structures by which (A, B)mutatesto
(A
,B
) and the corresponding function.
4. Analyze the relationship between the mutation mode and the alignment
mode of (A, B).
4
Super Pairwise Alignment
In Chap. 1, we introduced dynamic programming-based algorithms that are
used comprehensively in many fields. Next, we propose several modulus
structure-based and statistical decision-based algorithms. The key concept
giving rise to these algorithms was published in [90], and is referred to as
super pairwise alignment (SPA).
4.1 Principle of Statistical Decision-Based Algorithms
for Pairwise Sequences
In this section, we introduce pairwise alignment and the principle of SPA.
4.1.1 Uniform Alignment and Parameter Estimation
for Pairwise Sequences
The description of pairwise alignment was given in Chap. 1. The definitions
of minimum penalty alignment and uniform alignment were mentioned in
Sect. 1.2. The uniform alignment has been mentioned in Definition 23. In
this chapter, we focus on the solution for uniform alignment of pairwise se-
quences.
We continue to use the symbols given in (1.1). If sequence B is mutated
from A, and the mutation mode T is defined as in (3.55), then the solution for
uniform alignment of the pairwise sequence estimates the parameters {(i
k
,
k
),
k =1, 2, ··· ,k
a
} and k
a
.
If we estimate the parameter set T , which denotes the positions and lengths
of the mutations based on sequence (A, B), then sequence (A
,B
) constructed
by Definition 23 is an estimate for the uniform alignment of sequence (A, B).
Therefore, the key to solving the uniform alignment of the pairwise sequence is
110 4 Super Pairwise Alignment
knowing how to estimate the parameters in the mutation mode T basedonse-
quence (A, B). Then T in expression (3.55) is a group of statistical parameters
and
ˆ
T =
ˆ
i
k
,
ˆ
k
,k=1, 2, ··· ,
ˆ
k
a
(4.1)
is a set of statistics determined by (A, B),andanestimateoftheparameter
set T . The vital problem of uniform alignment of pairwise sequences is the
estimate of the parameters in T . The approach to solving this problem is
outlined as follows:
Sequential Estimate of Parameter Set T
The so-called sequential estimate for the parameter set T is briefly described
below:
1. To estimate the parameters in T alternately, we estimate (i
k
,
k
), k =
1, 2, ··· ,k
a
one after the other, that is, we estimate (i
k
,
k
) based on
(i
k
,
k
), k =1, 2, ··· ,k
− 1.
2. To estimate each (i
k
,
k
), we need not have the entire data of sequence
(A, B), but depend on only part of the data. Therefore, choosing the data
to use becomes one of the most important aspects of the statistical decision
algorithm.
3. The estimate of the parameter set T includes an estimate of the param-
eter k
a
.
Since this estimation procedure is identical to sequential estimation in statis-
tics, we call it the sequential estimation of the parameter set T .
The Locally Uniform Alignment of Sequence (A, B)
From the definition of uniform alignment in Definition 23, we generalize it to
formulate the definition of locally uniform alignment as follows:
Definition 24.
1. Let T be the shifting mutation mode given in expression (3.55), then
T
k
= {(i
k
,
k
),k=1, 2, ··· ,k
},k
≤ k
a
, (4.2)
is called the local shifting mutation of T .
2. Let B be the sequence mutated from A through shifting mutation mode T
and type-I and type-II mutations. If sequence B
k
is mutated from A
through shifting mutation T
k
and type-I and type-II mutations, then B
k
is called the local shifting mutation of A,andT
k
is corresponding local
shifting mutation mode.
3. If (A
,B
), (C
k
,D
k
) are two uniform alignmemt sequences of (A, B) and
(A, B
k
) respectively, then (C
k
,D
k
) is the locally uniform alignment se-
quence of (A
,B
).
4.1 Statistical Decision-Based Algorithms for Pairwise Sequences 111
The Outline of the Sequential Estimation for Parameter Set T
To estimate the parameter set T , we actually estimate (i
k
,
k
) and the locally
uniform alignment alternately. If the sequence pair (A, B) is given, then the
sequential estimation for parameter set T is summarized as follows:
1. Letting ¯a
0
=(a
1
,a
2
, ··· ,a
n
),
¯
b
0
=(b
1
,b
2
, ··· ,b
n
), we begin to estimate
(i
1
,
1
), and let (
ˆ
i
1
,
ˆ
1
) denote the estimation. The methods to select n
and to compute (
ˆ
i
1
,
ˆ
1
) will be detailed later.
2. After obtaining the estimation of the local shifting mutation T
k
,wedenote
it by
ˆ
T
k
similarly to (4.1), and then obtain the local alignment sequence
(
ˆ
C
k
,
ˆ
D
k
) based on
ˆ
T
k
and (A, B). Computation of the local alignment
sequence (
ˆ
C
k
,
ˆ
D
k
) will be described later.
3. Select vectors ¯c
t
=(c
t+1
,c
t+2
, ··· ,c
t+n
)and
¯
d
t
=(d
t+1
,d
t+2
, ··· ,d
t+n
)
from the local alignment sequence (
ˆ
C
k
,
ˆ
D
k
), so we can estimate
(
ˆ
i
k
+1
,
ˆ
k
+1
)of(i
k
+1
,
k
+1
). The method to select ¯c
t
and
¯
d
t
will be ex-
plained later in the block.
4.1.2 The Locally Uniform Alignment Resulting
from Local Mutation
Now, we detail the key steps 1–3 given in Sects. 4.1.1–4.1.3. We begin by
discussing the locally uniform alignment induced by a local mutation.
Let T
k
be a known local shifting mutation mode, which is defined by
expression (4.2). We can then obtain its local alignment sequence following
from the local shifting mutation mode T
k
discussed in Chap. 3. We recall
this process briefly, and comment on the specifics. For simplicity, we omit the
subscript k
of T
k
.
Some Symbols
The expansion or decomposition induced by a shifting mutation mode T are
addressed in Sect. 3.3.2. Here:
1. Let the decomposition of (C
,D
) expanded from (A, B) based on mode T
be
Δ
= {Δ
τ
,τ=0, 1, 2} =
δ
τ,k
,k=1, 2, ···,
ˆ
k
τ
,τ=0, 1, 2
. (4.3)
If we arrange the areas in Δ
in order, then we obtain Δ
= N
=
{1, 2, ··· ,n
}, and the decomposition expression is:
Δ
= {Δ
0
,Δ
1
,Δ
2
} =(δ
1
,δ
2
, ··· ,δ
2k
) , (4.4)
where k
is the number of times the shifting mutation occurred in T ,
n
= n
a
+ L
+
, L
+
=
k:
k
>0
k
. Then,
δ
2k−1
∈ Δ
0
,δ
2k
∈ Δ
1
∪ Δ
2
.
112 4 Super Pairwise Alignment
2. The decompositions of N
a
and N
b
are determined by the decomposition
of Δ
:
Δ
τ
= {Δ
τ,0
,Δ
τ,1
} =(δ
τ,1
,δ
τ,2
, ··· ,δ
τ,2k
τ
) ,τ= a, b , (4.5)
where k
a
,k
b
are the numbers of times the type-IV and type-III mutations
occurred in the shifting mutation mode T .Thus,k
a
+ k
b
= k
.
The intervals in (4.5) form a set produced by two steps: deletion Δ
1
,Δ
2
from Δ
, and then rearrangement of the rest of the intervals. If we denote
this set by
Δ
τ
= {i
τ,1
,i
τ,2
, ··· ,i
τ,2k
τ
,i
τ,2k
τ
+1
},τ= a, b , (4.6)
in which
i
τ,1
=0≤ i
τ,2
<i
τ,3
< ···<i
τ,2k
τ
≤ i
τ,2k
τ
+1
= n
τ
and
δ
τ,k
=[i
τ,k
+1,i
τ,k+1
]=(i
τ,k
+1,i
τ,k
+2, ··· ,i
τ,k+1
) ,
then by (3.78) and (3.82), we obtain the interrelationship between Δ
a
,Δ
b
,
and Δ
as follows:
i
a,2k−1
= j
0,2k−1
− L
1,k
,
i
a,2k
= j
0,2k
− j
0,2k−1
,
i
b,2k−1
= j
1,2k−1
− L
2,k
,
i
b,2k
= j
1,2k
− j
1,2k−1
,
(4.7)
where
L
τ,k
=
k
:δ
k
∈Δ
τ
,j
1,2k
−1
<j
0,2k−1
| δ
k
|,τ=1, 2 . (4.8)
3. Moreover, we have
C
= {a
Δ
a,0
,a
Δ
a,1
,b
Δ
b,1
},D
= {b
Δ
b,0
,b
Δ
b,1
,a
Δ
a,1
} , (4.9)
where the local vectors are arranged according to the order in Δ
,and
then
Δ
τ,0
= {δ
τ,1
,δ
τ,3
, ··· ,δ
τ,2k
−1
},
Δ
τ,1
= {δ
τ,2
,δ
τ,4
, ··· ,δ
τ,2k
},
where τ =1, 2andδ
τ,k
=[i
τ,k
+1,i
τ,k+1
].
Uniform Alignment Induced by the Shifting Mutation Mode T
From the shifting mutation mode T and sequence (A, B), we find the expanded
sequence (C
,D
), and we then obtain the alignment (A
,B
)of(A, B) through
the following steps:
1. Modify the sequence C
to become sequence A
, by replacing each com-
ponent in the subvector a
Δ
a,1
with virtual symbols.
4.1 Statistical Decision-Based Algorithms for Pairwise Sequences 113
2. Modify the sequence D
to become sequence B
, by replacing each com-
ponent in the subvector b
Δ
b,1
with virtual symbols.
Therefore, (A
,B
) is the uniform alignment if T =
ˆ
T
k
is the uniform esti-
mation for the local mutation mode T
k
,andthen(
ˆ
C
k
,
ˆ
D
k
), created by the
above process, is the uniform estimation statistic of (A, B) according to the
local mutation T
k
. We can estimate the next mutation position of (C
k
,D
k
)
based on (
ˆ
C
δ
2k
+1
,
ˆ
D
δ
2k
+1
).
4.1.3 The Estimations of Mutation Position and Length
We propose the statistical approach for estimating the mutation position and
length (i
k
,
k
) in this subsection.
Sliding Window Function of Sequence
The sliding window function of sequence A, B is defined by
w(A, B; i, j, n)=w
a
[i+1,i+n]
,b
[j+1,j+n]
=
1
n
n
k=1
w (a
i+k
,b
j+k
) , (4.10)
where w(a, b) is the Hamming matrix. Similarly, the sliding window function
of the stochastic sequence A
∗
,B
∗
is defined as
w(A
∗
,B
∗
; i, j, n)=w
a
∗
[i+1,i+n]
,b
∗
[j+1,j+n]
=
1
n
n
k=1
w
a
∗
i+k
,b
∗
j+k
. (4.11)
We will see later that the sliding window function is the same as the local
penalty function.
If B
∗
is mutated from A
∗
,wecanusew(A
∗
,B
∗
; i, j, n)toestimatethe
first mutation position of (A
∗
,B
∗
). Therefore, we can separate the random
variable pair (a
∗
i+k
,b
∗
j+k
) in expression (4.11) into two segments as follows:
The first segment: (a
∗
i+k
,b
∗
j+k
), k =1, 2, ···,n
1
is the segment without shifting
mutations, here b
∗
j+k
is the random variable mutated from a
∗
i+k
after type-I
and type-II mutations. The joint probability distribution is
P
r
a
∗
i+k
= b
∗
j+k
=1− , P
r
a
∗
i+k
= b
∗
j+k
= .
The second segment: (a
∗
i+k
,b
∗
j+k
), k = n
1
+1,n
1
+2, ··· ,n is the segment
with shifting mutation, here mutations III and IV occur in front of b
∗
j+k
and
a
∗
i+k
. The joint probability distribution is
P
r
a
∗
i+k
= b
∗
j+k
=
1
4
,P
r
a
∗
i+k
= b
∗
j+k
=
3
4
.
Let δ
1
=[1,n
1
],δ
2
=[n
1
+1,n]. Then we have
E{w(a
i+k
,b
j+k
)} =
⎧
⎨
⎩
, if (a
i+k
,b
j+k
) ∈ δ
1
,
3
4
, if (a
i+k
,b
j+k
) ∈ δ
2
.
(4.12)
114 4 Super Pairwise Alignment
If (A
∗
,B
∗
) are two independent stochastic sequences, then applying the
law of large numbers we have
w(A
∗
,B
∗
; i, j, n) ∼
n
1
n
+
3n
2
4n
, (4.13)
where n
τ
= ||δ
τ
|| is the length of the interval δ
τ
.
Statistics for the Sliding Window Function
From the computation formula for the sliding window function (4.13), we can
develop statistical estimates of (n
1
,n
2
). If the parameters , n are known,
following from (4.13) and n
1
+ n
2
= n we can estimate (n
1
,n
2
)as:
⎧
⎪
⎨
⎪
⎩
ˆn
1
= n
3 − 4w
3 − 4
= n
1 −
w −
0.75 −
,
ˆn
2
= n
w −
0.75 −
,
(4.14)
in which w = w(A
∗
,B
∗
; i, j, n) is calculated directly using (4.11) and (A
∗
,B
∗
).
Next,wediscusstheestimateof(n
1
,n
2
) when the parameter n is known
while is still to be determined. Let (A, B)beafixedsample.Wethenhave:
1. We assume that δ
1
and δ
2
are two fixed intervals, say,
δ
1
=[1, 2, ··· ,n
1
] ,δ
2
=[n
1
+1,n
1
+2, ··· ,n] ,
in which n
2
= n − n
1
. Then, a
i+δ
1
and b
j+δ
1
are segments without the
shifting mutations, and a
i+δ
2
and b
j+δ
2
are segments with shifting muta-
tions. Let
δ =(δ
1
,δ
2
)=(1, 2, ··· ,n) .
Then, (4.13) is also true and
w = w(A, B; i, j, n)=w(a
i+δ
,b
j+δ
) .
2. If we add a shifting parameter h to the penalty function, that is, if we
calculate
w
= w(A, B; i + h, j + h, n)=w(a
i+h+δ
,b
j+h+δ
) ,
in which 0 <h<nis a positive integer but not too large, then the
interval δ can be decomposed as:
δ =(δ
1
,δ
2
) = ((1, 2, ··· ,n
1
− h), (n
1
− h +1,n
1
− h +2, ··· ,n)) .
If a
i+h+δ
1
and b
j+h+δ
1
are segments without shifting mutations, and
a
i+h+δ
2
and b
j+h+δ
2
are segments with shifting mutations, the formula
for w
is found to be:
w
= w(A, B; i + h, j + h, n) ∼
(n
1
− h)
n
+
3(n
2
+ h)
4n
. (4.15)
4.2 Operation Steps of the SPA and Its Improvement 115
Therefore, following from (4.13) and (4.15) we obtain the following coupled
algebraic equations:
⎧
⎪
⎨
⎪
⎩
n
1
n
+
3n
2
4n
= w,
(n
1
− h)
n
+
3(n
2
+ h)
4n
= w
,
(4.16)
in which n
1
+ n
2
= n,andw, w
are constants that can be directly calcu-
lated based on A, B, i, j, n, h. Therefore, the value of n
1
, can be solved
using (4.16).
The simultaneous equations of (4.16) can be reduced to the following form:
4n
1
+3(n − n
1
)=4nw ,
4(n
1
− h) +3(n − n
1
+ h)=4nw
.
Thus, we have
⎧
⎪
⎨
⎪
⎩
=
3
4
+
n
h
(w − w
) ,
n
1
=
h
w
− w
3
4
− w
.
(4.17)
If we specify w, w
in (4.17) by w(A
∗
,B
∗
; i, j, n)andw(A
∗
,B
∗
; i + h,
j + h, n), respectively, then the results of (4.17) are denoted by ˆ, ˆn
1
,
which is the estimate of , n
1
.
4.2 Operation Steps of the SPA and Its Improvement
Continuing from the above section, we introduce the SPA for pairwise align-
ment. Let us begin by introducing the operation steps of SPA.
4.2.1 Operation Steps of the SPA
Let (A, B) be two fixed sequences. Every algorithm has the goal of estimating
all the parameters in the mutation mode T . The SPA is no exception. Specifi-
cally, we first select the important parameters n, h, θ, θ
,τ. Here, n is selected
according to the convergence of the law of large numbers or the central limit
theorem. Typically, we choose n =20, 50, 100, 150, etc. θ, θ
are selected based
on the error rate of the type-I and type-II mutations and the error rate of two
independently random variables. Thus, we choose 0 <θ<θ
< 0.75. For the
parameters h, τ as two local modifications, we choose them to be proportional
to n; typically, τ = αn, h = βn,0<α, β<0.5, etc. The SPA is described
below:
116 4 Super Pairwise Alignment
Step 4.2.1 Estimate the first mutation position i
1
in T :
1. Initialize i = j = 0 and calculate w(A, B; i, j, n). If
w(A, B; i, j, n)=w ≥ θ
,
then let
ˆ
i
1
= 0. This means the shifting mutation occurs at the
beginning of [1,n]. Otherwise, go to step 4.2.1-(2).
2. In Step 4.2.1, procedure 1, if w ≤ θ, meaning no shifting mutation
occurs in [1,n], we put the starting point forward and consider
i = j = n −τ. Next, we calculate the corresponding w(A, B; i, j, n).
If
w(A, B; i, j, n)=w ≤ θ,
then let i = j =2(n − τ) and repeat Step 4.2.1, procedure 2 until
w(A, B; i, j, n) >θ.Letk
1
be the integer satisfying the following
requirements:
w(A, B; i, j, n)=w ≤ θ
if i = j = k
1
(n − τ), and w(A, B; i, j, n)=w>θif i = j =
(k
1
+1)(n −τ ). Proceed to Step 4.2.1, procedure 3 or procedure 4.
3. For i = j =(k
1
+1)(n − τ), if w(A, B; i, j, n)=w ≥ θ
,thenset
ˆ
i
1
=(k
1
+1)(n − τ). Otherwise, go to Step 4.2.1, procedure 4.
4. Following Step 4.2.1, procedures 1–3, we have θ<w<θ
if i =
j =(k
1
+1)(n − τ). Therefore, for the same n, compute w
=
w(A, B; i + h, j + h, n). If w
>w,calculate
ˆ
i
1
according to (4.17).
Otherwise, repeat Step 4.2.1, procedures 1–4 for the larger h and n,
until w
>w.
Therefore, through the use of Step 4.2.1, we may estimate
ˆ
i
1
of i
1
.
Step 4.2.2 Estimate
1
based on the estimation
ˆ
i
1
of the first mutation po-
sition in T . Typically,
w
A, B;
ˆ
i
1
+ ,
ˆ
i
1
,n
,w
A, B;
ˆ
i
1
,
ˆ
i
1
+ , n
,=1, 2, 3, ··· .
If pair (
ˆ
i
1
+ ,
ˆ
i
1
)orpair(
ˆ
i
1
,
ˆ
i
1
+ )satisfiesw ≤ θ =0.3or0.4, where
w is its corresponding sliding window function, then this is the length
of the shifting mutation. Specifically:
1. If w(A, B;
ˆ
i
1
+ ,
ˆ
i
1
,n) ≤ θ,wenotethat
ˆ
1
= − and we insert
virtual symbols into sequence B following the position
ˆ
i
1
, while
keeping sequence A invariant.
2. If w(A, B;
ˆ
i
1
,
ˆ
i
1
+ , n) ≤ θ,wenotethat
ˆ
1
= and we insert
virtual symbols into sequence A following the position
ˆ
i
1
, while
keeping sequence B invariant.
Through the use of these two steps, we may estimate the local mutation
mode T
1
= {(i
1
,
1
)}, and its corresponding locally uniform alignment
(C
1
,D
1
). It is decomposed as follows:
C
1
=(C
1,1
,A
2,1
) ,D
1
=(D
1,1
,B
2,1
) .
4.2 Operation Steps of the SPA and Its Improvement 117
Denote the length of vector C
1,1
and D
1,1
by
ˆ
i
1
+ |
ˆ
1
|. Since there is no
shifting mutation occurring in the first n positions of A
2,1
,B
2,1
,welet
L
1
=
ˆ
i
1
+ |
ˆ
1
| + n be the starting point for next alignment.
Step 4.2.3 After obtaining the estimation (
ˆ
i
1
,
ˆ
1
), we continue to estimate i
2
basedon(C
1
,D
1
). We initialize i = j = L
1
and calculate w(A, B; i, j, n)
by repeating Step 4.2.1, procedures 1–4 to obtain the estimation
ˆ
i
2
for i
2
.
Step 4.2.4 Estimate
2
basedontheestimations
ˆ
i
1
,
ˆ
1
,
ˆ
i
2
. Here, we calculate
w
C
1
,D
1
;
ˆ
i
2
+ ,
ˆ
i
2
,n
,w
C
1
,D
1
;
ˆ
i
2
,
ˆ
i
2
+ , n
,=1, 2, 3, ··· .
We repeat Step 4.2.2 to get
ˆ
2
and the local alignment (C
2
,D
2
).
Step 4.2.5 Continuing the above process, we find the sequence (
ˆ
i
k
,
ˆ
k
)and
the corresponding sequence (C
k
,D
k
) for all k =1, 2, 3, ··· . The process
will terminate at some k
0
such that C
k
0
=(C
1,k
0
,A
2,k
0
)andD
k
0
=
(D
1,k
0
,B
2,k
0
) have shifting mutations occurring in (A
2,k
0
,B
2,k
0
). Let
L
k
0
denote the length of sequence C
1,k
0
,D
1,k
0
and i = j = L
k
0
.The
corresponding is the length of the shifting mutation if pair (
ˆ
i
k
0
+,
ˆ
i
k
0
)
or pair (
ˆ
i
k
0
,
ˆ
i
k
0
+ )satisfiesw ≤ θ,andthenw(C
k
0
,D
k
0
; i, j, n
) ≤ θ in
which n
is the shorter of the lengths of A
2,k
0
and B
2,k
0
.
Finally, we equalize the lengths of A
2,k
0
and B
2,k
0
.Inotherwords,if
the length of A
2,k
0
is shorter than that of B
2,k
0
, we insert several virtual
symbols at the end of A
2,k
0
so that its length is same as that of B
2,k
0
.
Example 17. The following two RNA sequences are drawn from 195 sRNAs
in [77]. We show how to align them using SPA:
(
E.co
: ugccuggcgg ccguagcgcg guggucccac cugaccccau gccgaacuca gaagugaaa
B.st
: ccuagugaca auagcggaga ggaaacaccc gucccauccc gaacacggaa guuaag
Here, n
a
=59,n
b
= 56. These two sequences seem disorderly and unsystem-
atic, but we can adapt them by inserting the virtual symbol “−” several times,
which is restricted so that the penalty functions of the two sequences are at
a minimum. By performing Steps 4.2.1–4.2.4, we obtain:
(
E.co
: ugccuggcgg ccguagcgcg guggucccac cugaccccau gccgaacuca gaagugaaa
B.st
: —ccuaguga caauagcgga gaggaaacac ccgucc-cau cccgaacacg gaaguuaag
The specific computational procedure is detailed as follows:
1. Let i = j =0,n = 15 and calculate the sliding window function of B.st
and E.co: w(A, B; i, j, n)=
14
15
=0.933 >θ
=0.6. Thus,
ˆ
i
1
=0.
2. For fixed
ˆ
i
1
=0,letn = 20 and calculate the sliding window functions
w(A, B;0, 1,n)=
13
20
=0.65 ,w(A, B;0, 2,n)=
18
20
=0.9 ,
w(A, B;0, 3,n)=
15
20
=0.75 ,w(A, B;0, 4,n)=
16
20
=0.8 ,
w(A, B;1, 0,n)=
16
20
=0.8 ,w(A, B;2, 0,n)=
7
20
=0.35 ,