跳到主要内容

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

Author

zephyr-zdz

Description

We have mm DNA sequences, and each sequence has the same length nn. The sequences contain four types of base: a, c, g, t. We wish to find all pairs of sequences where the last xx bases of one are identical to the first xx bases of the other (x<nx < n).

PP is an array of size mm, which holds one pointer to each sequence. It has already been sorted according to the lexicographic order of the sequences.

QQ is an array of size mm, which holds one pointer to each sequence. It has already been sorted according to the lexicographic order of the last xx bases of each sequence.

(1) Write pseudocode for an algorithm that outputs all pairs of sequences where the last xx bases of one are identical to the first xx bases of the other. It should exclude self-matches where one sequence has the same first and last xx bases. The running time of the algorithm should be linearly proportional to: mm plus the output size.

Your algorithm can use the following routines.

prefixSuffixCompare compares the length-xx prefix of its first argument to the length-xx suffix of its second argument, for example:

prefixSuffixCompare(P[i],Q[j],x)\mathrm{prefixSuffixCompare}(P[i], Q[j], x)

It returns 1-1 if the prefix is lexicographically less than the suffix, 00 if they are equal, and +1+1 if the prefix is greater. You can assume that this routine uses a constant amount of running time.

output outputs a pair of sequences, for example:

output(P[i],Q[j])\mathrm{output}(P[i], Q[j])

In these examples, 0i<m0 \le i < m and 0j<m0 \le j < m.

(2) Next, we wish to find all prefix-suffix matches of length x+1x + 1. Write pseudocode of an algorithm that fills in an array RR of size mm, so that it holds one pointer to each sequence, sorted according to the lexicographic order of the last x+1x + 1 bases of each sequence. The running time of the algorithm should be linearly proportional to mm.

Your algorithm can use this notation: Q[i][z]Q[i][z] is the base at position zz in the sequence pointed to by Q[i]Q[i], where 0z<n0 \le z < n.

题目描述

给定 mm 条长度均为 nn 的 DNA 序列,碱基取自 a, c, g, t 四种。希望找出所有满足「一条序列的后 xx 个碱基与另一条序列的前 xx 个碱基完全相同」的序列对(x<nx < n)。

已知两个大小为 mm 的指针数组:

  • PP:每条序列一个指针,已按序列整体的字典序排好序。
  • QQ:每条序列一个指针,已按每条序列末尾 xx 个碱基的字典序排好序。
  1. 写出输出所有上述序列对的算法伪代码。需排除「某条序列自身的前 xx 与后 xx 个碱基相同」而与自己配对的情形。算法运行时间须与 mm 加输出规模成线性关系。可使用两个子程序:prefixSuffixCompare(P[i], Q[j], x) 比较第一个参数长度为 xx 的前缀与第二个参数长度为 xx 的后缀,前缀较小返回 1-1、相等返回 00、较大返回 +1+1,可假定其运行时间为常数;output(P[i], Q[j]) 输出一个序列对。
  2. 接着希望求长度为 x+1x+1 的前后缀匹配。写出算法伪代码,把每条序列的指针填入大小为 mm 的数组 RR,使 RR 按每条序列末尾 x+1x+1 个碱基的字典序排列,运行时间须与 mm 成线性关系。可使用记号 Q[i][z]Q[i][z] 表示 Q[i]Q[i] 所指序列中位置 zz 处的碱基(0z<n0 \le z < n)。

Kai

Overview

The key observation is that the two arrays handed to us are precisely the two sorted orders the problem needs: PP groups the sequences by their first xx bases, and QQ groups them by their last xx bases. Part (1) is therefore not a search problem but a join problem — the same equality key appears, already sorted, on both sides — so a single linear merge pass locates every pair of matching blocks, and the only remaining work is emitting their cross product, whose cost is charged to the output size. This is exactly what the required running time of "mm plus the output size" is hinting at: it rules out testing all m2m^2 pairs, and it also tells us the answer may be far larger than mm, so the cross product must be written out directly rather than filtered afterwards. Part (2) extends the suffix key by one base while staying linear. Since the new key is "the base at position nx1n - x - 1" followed by "the old key", and QQ is already sorted by the old key, one stable counting sort on that single new base upgrades the order — one pass of an LSD radix sort. Stability is the whole point: it is what makes re-sorting from scratch unnecessary, and it is what keeps the cost at O(m)O(m) instead of O(mlogm)O(m \log m).

关键在于:题目给出的两个数组恰好就是所需的两种有序性——PP 按前 xx 个碱基把序列分组,QQ 按后 xx 个碱基把序列分组。因此 (1) 不是查找问题,而是连接(join)问题:同一个相等键已经在两侧分别排好序,于是一趟线性归并即可找出所有相匹配的块对,剩下的工作只是输出它们的笛卡尔积,其代价计入输出规模。题目要求的运行时间「mm 加输出规模」正是这个提示:它排除了枚举全部 m2m^2 对的做法,同时也说明答案规模可能远大于 mm,所以必须直接输出笛卡尔积,而不是先枚举再过滤。(2) 要在保持线性的前提下把后缀键延长一个碱基。新键等于「位置 nx1n - x - 1 处的碱基」接上「旧键」,而 QQ 已按旧键有序,故只需对这一个新增碱基做一次稳定的计数排序即可完成升级——也就是 LSD 基数排序的一趟。稳定性是关键:正是它使得无需从头重排,也正是它把代价保持在 O(m)O(m) 而非 O(mlogm)O(m \log m)

(1)

Idea

PP is sorted by the lexicographic order of the whole sequences, so sequences sharing the same first xx bases occupy a contiguous block of PP. Likewise, QQ is sorted by the lexicographic order of the last xx bases, so sequences sharing the same last xx bases occupy a contiguous block of QQ.

Hence it suffices to perform a merge join between the key "first xx bases" on PP and the key "last xx bases" on QQ: whenever the two keys agree, output the full cross product of the two blocks.

Pseudocode

i ← 0
j ← 0
while i < m and j < m:
r ← prefixSuffixCompare(P[i], Q[j], x)
if r = -1: // prefix of P[i] < suffix of Q[j]
i ← i + 1
else if r = +1: // prefix of P[i] > suffix of Q[j]
j ← j + 1
else: // match: locate the ends of both blocks
i2 ← i
while i2 < m and prefixSuffixCompare(P[i2], Q[j], x) = 0:
i2 ← i2 + 1
j2 ← j
while j2 < m and prefixSuffixCompare(P[i], Q[j2], x) = 0:
j2 ← j2 + 1
for a ← i to i2 - 1:
for b ← j to j2 - 1:
if P[a] ≠ Q[b]: // exclude a sequence paired with itself
output(P[a], Q[b])
i ← i2
j ← j2

Here output(P[a],Q[b])\mathrm{output}(P[a], Q[b]) denotes the pair in which the last xx bases of Q[b]Q[b] coincide with the first xx bases of P[a]P[a]. The test P[a]Q[b]P[a] \ne Q[b] compares pointers, i.e. it checks whether both entries refer to the same sequence.

Correctness

Let AsA_s be the set of sequences whose first xx bases equal ss, and BsB_s the set of sequences whose last xx bases equal ss. The answer is exactly

s{(p,q)pAs, qBs, pq}\bigcup_{s} \left\{ (p, q) \mid p \in A_s,\ q \in B_s,\ p \ne q \right\}

AsA_s appears as a contiguous block of PP and BsB_s as a contiguous block of QQ, and both keys are non-decreasing along their arrays. The merge join advances the side holding the smaller key, so it enumerates every pair of blocks sharing a common key ss, exactly once. For each such pair it scans the cross product and drops only the pairs of a sequence with itself, which matches the expression above.

Complexity

  • Every comparison advances ii or jj by one, so the branches other than else are taken at most 2m2m times in total.
  • Locating the block ends costs (i2i)+(j2j)(i_2 - i) + (j_2 - j) per matching block pair. The blocks are disjoint intervals of PP and of QQ, so this sums to at most 2m2m over all block pairs.
  • The double loop runs abab times for a block pair of sizes aa and bb. The only iterations that produce no output are those pairing a sequence with itself, of which there are at most min(a,b)\min(a, b) per block pair and at most mm in total. So the double loop runs at most (output size) + m+\ m times.

Altogether the running time is O(m+output size)O(m + \text{output size}), as required.

思路

PP 按序列整体的字典序排好序,因此前 xx 个碱基相同的序列在 PP 中占据一段连续区间(块);同理 QQ 按末尾 xx 个碱基的字典序排好序,因此后 xx 个碱基相同的序列在 QQ 中也占据一段连续区间。

于是只需对 PP 的「前 xx 碱基」这一键与 QQ 的「后 xx 碱基」这一键做归并连接(merge join):两键相等时,输出两个块的完整笛卡尔积。

上述伪代码中,output(P[a],Q[b])\mathrm{output}(P[a], Q[b]) 表示「Q[b]Q[b] 的后 xx 个碱基与 P[a]P[a] 的前 xx 个碱基相同」这一对;P[a]Q[b]P[a] \ne Q[b] 是指针比较,用于判断两项是否指向同一条序列。

正确性

设前 xx 碱基等于 ss 的序列集合为 AsA_s,后 xx 碱基等于 ss 的序列集合为 BsB_s,则答案恰为

s{(p,q)pAs, qBs, pq}\bigcup_{s} \left\{ (p, q) \mid p \in A_s,\ q \in B_s,\ p \ne q \right\}

AsA_sPP 的连续区间,BsB_sQQ 的连续区间,且两侧的键沿数组单调不减。归并连接总是推进键较小的一侧,因此不重不漏地枚举出所有共享同一键 ss 的区间对;对每个区间对遍历笛卡尔积,仅剔除同一序列自身配对的情形,与上式一致。

复杂度

  • 每次比较都使 iijj 前进 11,故 else 以外的分支合计至多执行 2m2m 次。
  • 定位块端点对每个匹配块对花费 (i2i)+(j2j)(i_2 - i) + (j_2 - j) 次;这些块分别是 PPQQ 上互不相交的区间,故所有块对合计至多 2m2m 次。
  • 对大小为 a×ba \times b 的块对,二重循环执行 abab 次;其中不产生输出的只有同一序列自身配对的情形,每个块对至多 min(a,b)\min(a, b) 个,全部块对合计至多 mm 个。故二重循环合计至多执行(输出规模)+ m+\ m 次。

综上,总运行时间为 O(m+输出规模)O(m + \text{输出规模}),满足题目要求。

(2)

Idea

The last x+1x+1 bases decompose into "the base at position z=nx1z = n - x - 1" followed by "the last xx bases". So the lexicographic order of the last x+1x+1 bases is the lexicographic order of the key pair

( Q[i][z], last x bases )(\ Q[i][z],\ \text{last } x \text{ bases}\ )

in which Q[i][z]Q[i][z] is the more significant key.

QQ is already sorted by the less significant key (the last xx bases), so one stable counting sort on the base at position zz suffices — this is exactly one step of an LSD radix sort. The alphabet has only the four bases a, c, g, t (a constant), so the counting sort runs in O(m)O(m) time. Note that x<nx < n guarantees z=nx10z = n - x - 1 \ge 0.

Pseudocode

z ← n - x - 1                    // 0-indexed position of the (x+1)-th base from the end

for each base β in {a, c, g, t}: // initialize the counters
count[β] ← 0

for i ← 0 to m - 1: // count the bases at position z
count[Q[i][z]] ← count[Q[i][z]] + 1

total ← 0 // prefix sums give each base its starting offset
for β in (a, c, g, t): // lexicographic order a < c < g < t
start[β] ← total
total ← total + count[β]

for i ← 0 to m - 1: // scanning Q in its own order keeps this stable
β ← Q[i][z]
R[start[β]] ← Q[i]
start[β] ← start[β] + 1

Correctness

Counting sort is stable: elements with equal keys keep their relative order from the input. The input QQ is ordered by the last xx bases, so entries sharing the same base at position zz retain that order inside RR. The starting offsets are computed by prefix sums taken in the lexicographic order a, c, g, t, so entries are ordered by the base at position zz. Therefore RR is sorted by (Q[i][z],last x bases)(Q[i][z], \text{last } x \text{ bases}), i.e. by the lexicographic order of the last x+1x+1 bases.

Complexity

Each loop performs O(1)O(1) iterations (four bases) or O(m)O(m) iterations, and each iteration takes constant time. The total running time is therefore O(m)O(m), as required.

Finally, running the algorithm of (1) on PP and this RR with x+1x+1 in place of xx yields all prefix-suffix matches of length x+1x+1.

思路

末尾 x+1x+1 个碱基可拆成「位置 z=nx1z = n - x - 1 处的碱基」加「末尾 xx 个碱基」。因此末尾 x+1x+1 碱基的字典序,等同于键对

( Q[i][z], 末尾 x 个碱基 )(\ Q[i][z],\ \text{末尾 } x \text{ 个碱基}\ )

的字典序,其中 Q[i][z]Q[i][z] 是更高位的键。

QQ 已按低位键(末尾 xx 碱基)排好序,故只需对位置 zz 的碱基做一次稳定的计数排序(counting sort)即可——这正是 LSD 基数排序的一步。字母表只有 a, c, g, t 四种(常数个),因此该计数排序为 O(m)O(m)。注意由 x<nx < n 可保证 z=nx10z = n - x - 1 \ge 0

正确性

计数排序是稳定排序:键相同的元素保持输入中的相对顺序。输入 QQ 按末尾 xx 碱基有序,故位置 zz 碱基相同的元素在 RR 中仍保持该顺序;而写入起始位置是按 a, c, g, t 的字典序做前缀和得到的,故元素按位置 zz 的碱基有序。于是 RR(Q[i][z],末尾 x 碱基)(Q[i][z], \text{末尾 } x \text{ 碱基}) 的字典序排列,即按末尾 x+1x+1 碱基的字典序排列。

复杂度

各循环分别执行 O(1)O(1) 次(四种碱基)或 O(m)O(m) 次,每次迭代为常数时间,故总运行时间为 O(m)O(m),满足题目要求。

最后,把这样得到的 RRPP 一起、以 x+1x+1 代替 xx 运行 (1) 的算法,即可求出所有长度为 x+1x+1 的前后缀匹配。