minimap and minimap2 原理介绍

编辑 / Bioinformatics / 发布于2022-08-02 / 更新于2023-03-16 / 阅读 404

Minimap and miniasm

符号定义

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

计算 minimizers

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

哈希函数: ϕ:ΣkZ\phi : \Sigma ^k \rightarrow \Zeta
ϕ(A)=0,ϕ(C)=1,ϕ(G)=2,ϕ(T)=3\phi(A)=0,\phi(C)=1,\phi(G)=2,\phi(T)=3,对于一个 k-mer s=a1...aks=a_1...a_k 定义
ϕ(s)=ϕ(a1)×4k1+ϕ(a2)×4k2+...+ϕ(ak)\phi(s) = \phi(a_1) \times 4^{k-1} + \phi(a_2) \times 4^{k-2} + ... + \phi(a_k)

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

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

h=ϕ(π(sik,r))=min{ϕ(π(sj+pk,r)):0p<w,r{0,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)M(s)(h,i,r)\in \Mu(s)(h,i,r)M(s)(h,i',r')\in \Mu(s'),就认为找到了一个 minimizer hit(h,x,i,i’) (匹配上了) 其中 x=rrx = r \oplus r' 。 h 是最小hash值,x表示相对链(同一条链值为0 或反向互补链为 1),i 和 i’ 是两个序列上的位置。

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

{(i1i1)(i2i2)<ϵ,if x=0(i1+i1)(i2+i2)<ϵ,if x=1 \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}

解释:

----------- i1i_1 -----i2i_2 —>

i1i'_1 ------i2i'_2 ----------------->

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

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

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

Assembly graph

给了几个定义

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

Minimap2

Materials and methods

minimap2 遵循 seed-chain-align 的流程

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

chaining

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

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

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

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

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

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

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

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

Backtracking