6 Minseed Algorithm

下面把论文 第 6 节 MINSEED ALGORITHM 拆开讲清楚。
这一节讲的是 SeGraM 中负责 seeding（种子查找） 的算法：MinSeed。
它的目标是：为每条 query read 快速找到候选映射位置（子图），交给后面的 BitAlign 做精细比对。

---

一、核心思想：为什么用 minimizer？

MinSeed 基于 Minimap2 的 minimizer-based seeding 算法（mm_sketch）。
要理解它，先要明白什么是 minimizer。

1. 什么是 minimizer？

论文定义：

A $< w,k >$ minimizer is the smallest $k$-mer in a window of $w$ consecutive k-mers (according to a scoring mechanism), for subsequences of length $k$.

翻译：
一个 $< w,k >$ minimizer 是长度为 $w$ 的连续 k-mer 窗口中，按某种评分机制选出的最小的那个 k-mer。

• k：k-mer 的长度，比如 k=3；
• w：窗口大小，表示窗口里包含多少个连续的 k-mer；
• 窗口：在序列上滑动，每次包含 $w$ 个相邻的 k-mer；
• 最小：按字典序或其他排序规则，选最小的那个。

2. 为什么用 minimizer 而不是所有 k-mer？

论文说：

Compared to using the full set of k-mers, using only the $< w,k >$ minimizers decreases the storage requirements of the index (by a factor of $2/(w+1)$) and speeds up index queries.

翻译：
与使用全部 k-mer 相比，只用 minimizer 可以把索引存储量减少约 $2/(w+1)$ 倍，并加速索引查询。

通俗理解：
如果把 read 切成所有 k-mer，数量很多，索引会很大。
但每个窗口只选一个最小的 k-mer，就能大幅减少需要索引和查询的 k-mer 数量。
同时，minimizer 保证：如果两条序列有足够长的精确匹配，它们至少会共享一个 minimizer。
所以用 minimizer 做种子，既省空间又不会漏掉真正的匹配。

3. 一个简单例子

论文 Figure 8 给了例子：

序列：A G T A G C A ...
窗口 w=5，k=3

第一个窗口包含 5 个连续的 3-mer：

位置	k-mer
1	A G T
2	G T A
3	T A G
4	A G C
5	G C A

按字典序排序，最小的 k-mer 是 AGC。
所以第一个窗口的 minimizer 就是 AGC。

然后窗口向右滑动一个位置，重复这个过程，选出每个窗口的 minimizer。

---

二、MinSeed 算法流程

论文说：

The MinSeed algorithm starts by computing the minimizers of a given query read.

翻译：
MinSeed 算法首先计算给定 query read 的 minimizer。

然后按顺序做以下几件事：

步骤 1：计算 query read 的 minimizer

• 对 query read 滑动窗口，每个窗口选一个 minimizer；
• 把所有 minimizer 收集起来。

优化：
论文提到：

While a naive way to compute the minimizers is to use a nested loop ... we can eliminate the inner loop by caching the previous minimum k-mers within the current window. The single-loop algorithm has a complexity of $O(m)$, where $m$ is the length of the query read.

翻译：
朴素方法用嵌套循环：外层遍历窗口，内层找窗口内最小 k-mer。
但可以缓存上一个窗口的最小 k-mer，消除内层循环，使复杂度降为 $O(m)$，其中 $m$ 是 read 长度。

通俗理解：
窗口滑动时，很多 k-mer 是重复的，不需要每次都重新比较。
用一个高效的单循环就能在线性时间内找出所有 minimizer。

步骤 2：查询每个 minimizer 的出现频率

• 用 minimizer 作为键，去主存的哈希表索引里查它的出现频率（即它在参考图中出现了多少次）。

步骤 3：过滤高频 minimizer

论文说：

A minimizer is discarded if its occurrence frequency in the reference genome is above a user-defined threshold (pre-computed for each chromosome in order to discard the top $0.02\%$ most frequent minimizers), in order to reduce the number of seed locations that are sent to the alignment step.

翻译：
如果某个 minimizer 在参考基因组中的出现频率超过用户定义的阈值，就丢弃它。
阈值是预先为每条染色体计算的，目的是丢弃最频繁的 top 0.02% 的 minimizer，从而减少送到 alignment 步骤的 seed 数量。

为什么？
高频 minimizer 通常出现在重复区域，会产生大量候选位置。
如果每个都做 alignment，计算量会爆炸。
所以先丢掉这些“烂大街”的 minimizer。

步骤 4：获取剩余 minimizer 的 seed 位置

• 对保留下来的 minimizer，从索引里取出它们在参考图中的所有精确匹配位置（即 seed locations）；
• 每个 seed location 包含：节点 ID + 节点内偏移。

步骤 5：计算候选子图（左右位置）

论文说：

After fetching all seed locations ... MinSeed calculates the leftmost and rightmost positions of each seed, using the node ID and relative offset of the seed location along with the relative offset of the corresponding minimizer within the query read.

翻译：
取完所有 seed 位置后，MinSeed 计算每个 seed 的最左位置和最右位置，用到：

• seed 位置（节点 ID + 节点内偏移）；
• 对应 minimizer 在 query read 中的相对偏移；
• 错误率 E。

Figure 9 展示了计算方法：

query read:    0 ... a ... b ... m-1
graph ref:     x ... c ... d ... y

• a：minimizer 在 read 中的起始位置；
• b：minimizer 在 read 中的结束位置；
• c：seed 在参考中的起始位置；
• d：seed 在参考中的结束位置；
• m：read 长度；
• E：错误率。

左扩展 = $a \times (1+E)$
右扩展 = $(m - b - 1) \times (1+E)$

所以：

• 候选区域左边界：$x = c - a \times (1+E)$
• 候选区域右边界：$y = d + (m - b - 1) \times (1+E)$

通俗理解：
minimizer 只是 read 的一小段。
read 的其余部分可能落在 minimizer 左边或右边。
考虑到测序错误和变异，左右各留出一定余量（用错误率 E 放大），就得到整个 read 可能匹配的参考区域，也就是候选子图。

步骤 6：把候选子图交给 BitAlign

• 从主存中取出这些候选子图涉及的图节点；
• 写入 BitAlign 的 input scratchpad，等待精细比对。

---

三、MinSeed 的设计要点

1. 单循环复杂度 O(m)

计算 minimizer 用单循环，线性时间，很快。
因为 seeding 需要为每条 read 反复做，所以效率很重要。

2. 频率过滤代替复杂过滤

论文在 11.4 节提到：

MinSeed does not implement a filtering mechanism. ... MinSeed does not decrease the sensitivity ... since both MinSeed and the baseline software tools implement the same optimization of discarding the seeds that have higher frequency than the threshold.

翻译：
MinSeed 不实现额外的过滤机制，只做频率过滤。
它不会降低整体灵敏度，因为基线的软件工具也做了同样的高频种子丢弃优化。

意思：
MinSeed 不做 chaining、clustering 等复杂过滤，但通过丢弃高频 minimizer，已经能大幅减少候选数量。
剩下的候选交给 BitAlign，因为 BitAlign 很快，所以整体性能仍然很好。

3. 不降低灵敏度

灵敏度 = 保留真正正确 seed 的能力。
MinSeed 只丢弃高频 minimizer，这些通常来自重复区域，丢弃它们不会漏掉真正的独特匹配。
所以灵敏度与基线工具相当。

---

四、MinSeed 与 BitAlign 的衔接

MinSeed 的输出是候选子图，不是最终比对结果。
这些子图被送到 BitAlign，由 BitAlign 做近似字符串匹配和编辑距离计算。
因为 MinSeed 已经把候选范围缩小了很多，BitAlign 只需要处理少量子图，整体速度就上去了。

论文在 11.4 节举了例子：

• 长读数据集：GraphAligner 用过滤/链接把 seed 从 77M 降到 48k；MinSeed 只降到 35M，但 SeGraM 仍然比 GraphAligner 快，因为 BitAlign 效率极高。
• 短读数据集：GraphAligner 从 828k 降到 11k；MinSeed 只降到 375k，SeGraM 仍显著胜出。

这说明：MinSeed 不需要做很激进的过滤，因为 BitAlign 足够快，能消化更多候选。

---

五、总结：第 6 节的核心逻辑

1. MinSeed 基于 minimizer：每个窗口选最小 k-mer，减少索引和查询量。
2. 算法步骤：
   • 计算 read 的 minimizer（O(m) 单循环）；
   • 查频率，丢弃高频 minimizer；
   • 取剩余 minimizer 的 seed 位置；
   • 根据 minimizer 在 read 中的位置和错误率，计算候选子图的左右边界。
3. 频率过滤：丢弃 top 0.02% 最频繁的 minimizer，减少候选数量。
4. 不做复杂过滤：因为 BitAlign 很快，不需要 MinSeed 做过度过滤。
5. 输出候选子图：交给 BitAlign 做精细比对。

一句话：
MinSeed 用 minimizer 快速找出 read 在参考图中的候选位置，通过频率过滤去掉重复区域的种子，再根据错误率扩展出候选子图，最后把这些子图交给 BitAlign。它用线性时间完成 seeding，为后续 alignment 大幅缩小了搜索范围。