原文链接:Reducing storage requirements for biological sequence comparison | Bioinformatics | Oxford Academic (oup.com)
论文封面
对于序列比对问题,可以细分为两种,
- 一种是叫做align(比对),如bsalign、ksw等
- 一种是叫做map(映射),如bwa、minimap2、bitmapper等
对于align问题,往往需要知道两条或者多条序列比对的分数,来借此衡量序列的相似程度。对于map问题,则是将短的序列映射到长的序列上,一个很常见的应用是在有参组装基因组时,将短的reads映射到参考基因组上,从而将reads定位在基因组的相对位置上。align对于精度要求高,而map则对于精度要求没有那么高。
如何做map!
将150bp的reads map到基因组上
一个较为精确的方法是通过align的 方法,计算一个打分矩阵,通过回溯来确定reads在基因组(ref)上的位置。这种方法可以完成比对的目标,每map一条reads,需要计算一个巨大的打分矩阵,而实际上大部分的计算是没有必要的,150bp的长度在基因组水平上简直是沧海一粟,所以这种使用纯align的方法可行,但实际上不可行。原因有二,一是不需要太高的精度(当然精度高也可以),二是无用计算占比大,map速度慢(主要原因)。
seed and extend
既然上面的方法在实现上不同,主要原因是计算量大,且大部分是无效计算。那么有没有办法减少无效计算哪?
基于seed-and-extend的方法可以做到这一点,有两条序列1和2,外面首先在两个序列中各找到一个匹配的seed(图中绿色的部分),然后将两条序列的seed对齐,在以seed为中心,向两侧延伸进行align,当align的分数小于一定值,或者是两侧扩展超过一定程度时就停止,这样就可以做到即减少计算量,同时也保证了一定的精度。
何为seed
上面的seed可以被叫做k-mer,一种在序列上进行平移的短序列,这里的k为短序列的长度,比如,下面的4-mer就是从序列的第一个碱基开始,以4的窗口长度进行向后平移,每次位移的长度为1,在这个过程中便产生了一系列的4-mer。
k-mer
使用kmer可以将一条序列转化为若干个种子,从而将不同的序列进行粗对齐(其中需要考虑共线性)。但是kmer存在一个问题,就是占用空间太大,比如,上面图片中的序列长度为10,产生的4-mer为7个,也就是有28个碱基。当序列的长度为n时,产生的k-mer数位n-k+1,所有的kmer加起来有k*(n-k+1),相比于原序列空间几乎增长了k倍。
minimizer
通过在一定的窗口大小内,有指向的选择1个kmer,这一个kmer被称为minimizer。如下面的例子,我们规定 A>C>G>T,那么我们选择最大的一条kmer,即为绿色的AGCTGAGCT,通过这种人为的规定什么样的kmer才是大,什么样的kmer才是小,我们可以在相同的序列中得到相同的minimizer。
minimizer的选择依赖于两个参数
- k :即kmer的长度
- w:即w个相邻的kmer中选择一个kmer作为minimizer
不同的参数造成的结果不同
- k=4, w=3
k=4, w=3
- k=4, w=4
k=4, w=4
- k=4, w=5
k=4, w=5
- k=4, w=6
k=4, w=6
通过选择不同的k和w,可以发现,当w>k时,会造成选择出来的minimizer没有覆盖到序列中中间部分,当w<=k时,可以保证选择出的minimizer一定可以覆盖到中间部分,不会造成中间位置某个地方产生空缺。所以在实际的应用中推荐w<=k。仔细观察,还可以发现,无论k和w是什么关系,都有可能造成序列的开头和结尾部分产生空缺,所以一个有效的方法是对开头和结尾做特殊处理。
有指向的选择1个kmer
在上一段中我们提到了minimizer为在w个相邻的kmer中“有指向的选择1个kmer”,这里的有指向值得玩味。我们可以将“有指向抽象为一个数学函数”,比如F(kmer_1, kmer_2, ..., kmer_w)
,只要是w个相邻的kmer所对应的序列是相同的的那么F(kmer_1, kmer_2, ..., kmer_w)
得到的minimizer也一定是一样的。这样做可以保证,相同的一段序列总能得到相同的minimizer,这对我们最开始说使用seed将两条序列做对齐十分有帮助。但是无论是什么样的F,都只能从w个kmer中选择一个minimizer,如果仅仅依靠minimizer相同就推断序列相同是错误的,还需要依靠其他的信息,比如minimizer的共线性信息(minimap中的技术)等。