Skip to content

東京大学 新領域創成科学研究科 メディカル情報生命専攻 2024年1月実施 問題12

Author

zephyr

Description

ナノポアシークエンサーは DNA 配列を電流値の変化を通して読み取ることができる。 ある DNA 配列を読み取ったところ、時間 \(t = 1, 2, \ldots , n\) でそれぞれ読み取られた電流値列は \(\{A_t\} \ (0 \leq A_t \leq 1)\) であった。 この DNA 配列には基準配列が含まれていることが分かっており、基準配列を単独で読み取ったときの電流値列は時間 \(t = 1, 2, \ldots, m\) に対して \(\{B_t\} \ (0 \leq B_t \leq 1)\) であった（基準波形）。 このとき、\(A_t\) の中で基準配列に相当する時間区間を見つけたい。

DNA 配列を読み取る速度にはばらつきがあり急な変動があるため、基準配列を探す場合には時間軸方向へのずれを許したい。 このため、\(1 \leq p_1 \leq p_2 \leq \cdots \leq p_m \leq n\) (ただし \(1 leq t < m\) に対して \(p_{t+1} - p_t \leq W\), \(W\) は正の整数) が存在して \(B_t\) と \(A_{p_t}\) が対応していると仮定する。 読み取りノイズ等が無ければ \(B_t = A_{p_t}\) が成立して欲しいが、実際には読み取りノイズや電流値読み取りタイミングのずれを考慮して、相違度 \(\sum_{t=1}^m (B_t - A_{p_t})^2\) が最小値をとるような対応関係 \(p_1, p_2, \ldots, p_m\) が基準波形とナノポアシークエンサー出力の一部分との正しい対応関係を示していると仮定する。

(1) \(m = 1\) とする。相違度が最小となる \(p_1\) を１つ求めるアルゴリズムを示せ。

(2) \(m = 2\) とする。相違度が最小となる \(p_1\), \(p_2\) の組を１つ求めるアルゴリズムを示せ。

(3) \(W\) はどのような物理条件に対応していると考えられるか。１行で説明せよ。

(4) \(D_{i,j}\) は \(B_i\) と \(A_j\) が対応しているときの、基準配列 \(\{B_1, B_2, \ldots, B_i\}\) に対する \(\{A_1, A_2, \ldots, A_j\}\) の相違度の最小値である。\(D_{x,y} \ (x < i, y \leq j)\) を用いて \(D_{i,j}\) を表せ。

(5) 基準配列 \(\{B_t\}\) に対する \(\{A_t\}\) の相違度の最小値を計算するアルゴリズムを示し、その時間計算量を \(n\), \(m\), \(W\) を用いて示せ。

Kai

1. Finding \(p_1\) for \(m = 1\)

When \(m = 1\), the problem reduces to finding a single time point \(p_1\) such that the dissimilarity \(\sum_{t=1}^1 (B_t - A_{p_t})^2\) is minimized. The algorithm can be described as follows:

  1. Initialize \(p_1\) to 1.
  2. For each \(i = 1, 2, \dots, n\), compute the dissimilarity \(d_i = (B_1 - A_i)^2\).
  3. Find the index \(i\) that minimizes \(d_i\).

The resulting \(p_1\) will be the index that minimizes the dissimilarity.

Algorithm

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
Input: A sequence {A_t}, reference value B_1
Output: Index p_1 that minimizes dissimilarity

1. Initialize min_dissimilarity = infinity
2. Initialize p_1 = 1
3. For i = 1 to n do
     d_i = (B_1 - A_i)^2
     If d_i < min_dissimilarity then
         min_dissimilarity = d_i
         p_1 = i
4. Return p_1

2. Finding \(p_1, p_2\) for \(m = 2\)

When \(m = 2\), the problem involves finding two indices \(p_1\) and \(p_2\) such that \(1 \leq p_1 \leq p_2 \leq n\) and \(p_{t+1} - p_t \leq W\). The goal is to minimize the dissimilarity \(\sum_{t=1}^2 (B_t - A_{p_t})^2\). The algorithm can be described as follows:

  1. Initialize minimum dissimilarity to infinity and \(p_1, p_2\) to 1.
  2. For each \(i = 1, 2, \dots, n\), do:
  3. For each \(j = i, i+1, \dots, \min(i + W, n)\), do:
    • Compute the dissimilarity \(d_{i,j} = (B_1 - A_i)^2 + (B_2 - A_j)^2\).
    • If \(d_{i,j} < \text{min\_dissimilarity}\), update \(p_1 = i\), \(p_2 = j\), and \(\text{min\_dissimilarity} = d_{i,j}\).

Algorithm

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
Input: A sequence {A_t}, reference values {B_1, B_2}, window W
Output: Indices p_1, p_2 that minimize dissimilarity

1. Initialize min_dissimilarity = infinity
2. Initialize p_1, p_2 = 1
3. For i = 1 to n do
     For j = i to min(i + W, n) do
         d_ij = (B_1 - A_i)^2 + (B_2 - A_j)^2
         If d_ij < min_dissimilarity then
             min_dissimilarity = d_ij
             p_1 = i
             p_2 = j
4. Return p_1, p_2

3. Physical Meaning of \(W\)

The parameter \(W\) represents the maximum allowable fluctuation in the time axis when matching the reference signal to the nanopore sequencer output. It accounts for variations in the speed at which the DNA is read, allowing for some flexibility in aligning the sequences.

4. Recursive Formula for \(D_{i,j}\)

The minimum dissimilarity \(D_{i,j}\) can be computed using the following recurrence relation:

\[ \begin{aligned} D_{i,j} &= \min_{x < i, y \leq j} \{D_{x,y} + (B_i - A_j)^2 \} \quad \text{subject to } j - y \leq W \\ &= \min_{\max(1, j-W) \leq k < j} {D_{i-1,k} + (B_i - A_j)^2} \end{aligned} \]

5. Algorithm for Minimum Dissimilarity for \(m > 2\)

To calculate the minimum dissimilarity between the sequence \(\{A_t\}\) and the reference signal \(\{B_t\}\), we use a dynamic programming approach.

  1. Initialize a 2D array \(D\) of size \(m \times n\) with infinity, where \(D_{i,j}\) represents the minimum dissimilarity for the first \(i\) elements of the reference sequence mapped to the first \(j\) elements of the output sequence.
  2. Set \(D_{1,j} = (B_1 - A_j)^2\) for all \(j = 1, 2, \dots, n\).
  3. For each \(i = 2, 3, \dots, m\), do:
  4. For each \(j = 1, 2, \dots, n\), do:
    • Set \(D_{i,j} = \min_{x = \max(1, j - W)}^{j-1} \{D_{i-1,x} + (B_i - A_j)^2\}\).

Time Complexity

The time complexity of this algorithm is \(O(m \cdot n \cdot W)\), where:

  • \(m\) is the length of the reference signal,
  • \(n\) is the length of the sequence \(\{A_t\}\),
  • \(W\) is the maximum allowed shift between adjacent points.

Knowledge

难点思路

本题的主要难点在于第 4 和第 5 问，需要正确地推导出动态规划的递推关系，并设计出高效的算法。关键是要理解如何处理时间轴偏移的约束条件，以及如何在保证正确性的同时优化算法的时间复杂度。

解题技巧

对于序列对齐问题，动态规划是一种有效的解决方案。特别是在处理时间序列中的波动和对齐问题时，可以灵活调整对齐窗口 \(W\) 来达到最优匹配。

重点词汇

  • Dissimilarity: 相异度
  • Reference Signal: 参考信号
  • Window: 窗口

参考资料

  1. Dynamic Time Warping Algorithm: Introduction to the dynamic time warping algorithm, Wikipedia.
  2. Time Series Analysis and Its Applications, Third Edition, Robert H. Shumway and David S. Stoffer.