Minimap and miniasm

符号定义

$\sum = {A,C,G,T}$ 为核苷酸的字符表。对于 $a \in \sum$ ， $\overline a$ 是 $a$ 的的互补字符（A-T，C-G）。字符串 $s=a_1a_2...a_n$ 也叫做DNA序列 $|s|=n$ 是它的长度，它的反向互补 $\overline s = \overline{a_1a_2...a_n} = \overline a_n \overline a_{n-1} ... \overline a_1$ 。定义链函数 $\pi: \sum ^* \times \{0,1\} -> \sum^*$ ，就像 $\pi(s,0) = s$ 和 $\pi(s,1) = \overline s$ 。把长度为k的DNA序列成为k-mer，用 $s_i^k = a_i...a_{i+k-1}$ 标记一个从 i 开始的，长度为 k 的，s 的子序列。 $\sum ^k$ 是所有 k-mer 的集合。

计算 minimizers

一个字符串的 (w,k)-minimizer 是在 w 个连续的 k-mer 周围窗口中最小的那个 k-mer。

哈希函数： $\phi : \Sigma ^k \rightarrow \Zeta$ 。
$\phi(A)=0,\phi(C)=1,\phi(G)=2,\phi(T)=3$ ，对于一个 k-mer $s=a_1...a_k$ 定义
$\phi(s) = \phi(a_1) \times 4^{k-1} + \phi(a_2) \times 4^{k-2} + ... + \phi(a_k)$

双链 (w,k, $\phi$ ) minimizer $|s| \ge w+k-1$ ，可以写作 (h,i,r), 其中，h是 minimizer 的值，i是处于序列上的位置(这个k-mer在序列上的开始位置)，r是位于那一条链上（反向互补链）

$max(1,i-w+1) \le j \le min(i,|s| -w -k +1)$

$h=\phi(\pi(s_i^k,r))=min\{\phi(\pi(s_{j+p}^k,r')):0\le p < w, r' \in\{0,1\}\}$

$\Mu $ 作为字符串 s 的 minimizer 集合。

Indexing

把所有目标序列的 minimizers 放到 hash table 中，minimizer 的 hash 值作为 key， value 是目标序列索引，minimizer的位置和链的方向的集合。

Mapping

两个序列 s 和 s’, 如果存在 $(h,i,r)\in \Mu(s)$ 和 $(h,i',r')\in \Mu(s')$ ，就认为找到了一个 minimizer hit(h,x,i,i’) (匹配上了) 其中 $x = r \oplus r'$ 。 h 是最小hash值，x表示相对链（同一条链值为0 或反向互补链为 1），i 和 i’ 是两个序列上的位置。

另外如果两个 minimizer hits $hit(h_1,x,i_1,i'_1)$ 和 $hit(h_2,x,i_2,i'_2)$ 能满足

\begin{cases} |(i_1 - i'_1) -(i_2-i'_2)| <\epsilon , & \text{if x=0} \\ |(i_1 + i'_1) -(i_2 + i'_2)| <\epsilon , & \text{if x=1} \\ \end{cases}

解释：

----------- $i_1$ ----- $i_2$ —>

— $i'_1$ ------ $i'_2$ ----------------->

如上所示，上述的计算是为了判断两个 hits 分别在两条序列上距离的差值，如果两个hits的差值小于 $\epsilon$ 就称为 $\epsilon-away$

则称他们是 $\epsilon-away$ 的。 $\epsilon-away$ 的 hits 在宽度 $\epsilon$ 内近似的共线性。
给定一个 $hit(h_1,x,i_1,i'_1)$ 集合，可以通过聚类(对于 x=0 i-i’，对于 x=1 i+i’ )来确定长的共线性匹配

这个算法里不同 t 值的 hit 不会被分到同一个集合中，感觉这个地方有点问题。

Assembly graph

给了几个定义

w contains v : 如果 v 能够比对到 w 的子序列里
w overlaps v ：如果 v 的后缀序列（从内部某个位置开始到结束的子序列）和 w 的前缀序列可以互相比对。写作 $v->w$

Minimap2

Materials and methods

minimap2 遵循 seed-chain-align 的流程

首先收集参考序列上的 minimizers，然后把它放进 hash table 里，minimizer 的 hash 值作为 key，value 是这个 minimizer值所在的参考序列的位置列表（多个位置的 minimizer 的值可能是相同的。
然后对于每个查询序列，把多有查询序列的 minimizers 作为种子，找到在参考序列里面相匹配的（anchors），然后找到所有共线性的 anchors 的集合作为 chains，
如果需要单个碱基级别的比对信息，minimap2 使用动态规划从 chains 的端点延伸并把 chains 中的相邻 anchors 之间的区域连接起来。

chaining

anchor ： $(x,y,w)$ 表明参考序列上的 $[x-w+1,x]$ 和查询序列（query read）上的 $[y-w+1,y]$ 区域匹配（minimizer值相同）。

把所有 anchors 按照在参考序列的结束位置 x 排序， $f(i)$ 作为列表里第 i 个 anchor 最大的 chaining score。

f(i)=max\{max \ i>j≥1\ \{f(j)+α(j,i)−β(j,i)\},w_i\}

其中 $α(j,i)=min\{min\{y_i−y_j,x_i−x_j\},w_i\}$ ，是两个 anchors 之间匹配的碱基数量（感觉不管有没有匹配都加上了）。 $β(j,i)>0$ 是 gap cost （减分）。如果 $y_j≥y_i\ or\ max\{y_i−y_j,x_i−x_j\}>G$ 它为 $\infty$ ，否则

β(j,i)=γ_c((y_i−y_j)−(x_i−x_j))

这里的 $(y_i−y_j)−(x_i−x_j)$ 就和 Mapping 那里是一样的。

y_c = \begin{cases} 0.01\times \overline w \times |l| + 0.5log2 & (l=0)\\ 0 & (l=0)\\ \end{cases}

其中 $\overline w$ 是平均种子长度

minimap and minimap2 原理介绍