BWT算法浅析-(一)
一、前言
BWT算法对数据量大的数据来说,压缩比较有效,目前已经广泛用于生物学数据的存储中,例如bowtie就使用到了这一策略。虽然之前在大学期间学过该算法,但是有点遗忘了,故重新复习一番。
二、什么是BWT?
严格来说,BWT并不是压缩算法,它仅是把需要压缩的字符重新排列,将相同的字符排在一起或者相邻位置,使其更易被压缩。
三、BWT原理
BWT算法主要分为两个部分,编码以及解码。
3.1、编码
编码过程主要分为3个步骤:
1. 循环待压缩字符,构建字符矩阵M; 2. 按照字典序重排矩阵M; 3. 提取重排后的矩阵M的最后一列L;
下面我们举个例子来说明,假如我们需要压缩字符串S是ACGTAA。
首先,我们在该字符后添加一个字符$(该字符不存在于字符串S中,且其ASCII码小于字符串S中的任意一个字符,思考一下为什么呢?),然后将当前字符串的第一个字符移到最后一位,形成一个新的字符串,再将新的字符串的第一位移到最后一位形成另一个新的字符串,就这样不断循环这个过程,直到字符串循环完毕(即$处于第一位),这样就形成了一个基于原字符串的字符矩阵。这里我们使用python代码来进行实现:
def rotations(s):
s_mat = []
for i in range(0, len(s)):
tmp = s[i:]+s[:i]
s_mat.append(tmp)
return s_mat
接下来就是按照字典序排序矩阵M,这个比较简单,就不多说了,在python中直接sort即可。矩阵如下,其中F列为矩阵第一列,L列为矩阵最后一列
ID | rotations | sort | F | L |
---|---|---|---|---|
1 | ACGTAA$ | $ACGTAA | $ | A |
2 | CGTAA$A | A$ACGTA | A | A |
3 | GTAA$AC | AA$ACGT | A | T |
4 | TAA$ACG | ACGTAA$ | A | $ |
5 | AA$ACGT | CGTAA$A | C | A |
6 | A$ACGTA | GTAA$AC | G | C |
7 | $ACGTAA | TAA$ACG | T | G |
3.2、解码
解码主要是通过L列还原到原来字符串,其主要原理是通过LF列对应关系进行还原。为了方便说明,我们给每个出现的字符一个编号,表示当前字符出现了几次,原字符串为A0C0G0T0A1A2。
Sort | F | L |
---|---|---|
$0A0C0G0T0A1A2 | $0 | A2 |
A2$0A0C0G0T0A1 | A2 | A1 |
A1A2$0A0C0G0T0 | A1 | T0 |
A0C0G0T0A1A2$0 | A0 | $0 |
C0G0T0A1A2$0A0 | C0 | A0 |
G0T0A1A2$0A0C0 | G0 | C0 |
T0A1A2$0A0C0G0 | T0 | G0 |
从上表中可以看出,F列中对应字符出现的顺序与L列中一致,例如:F列中A2A1A0与L列中A的出现顺序一致。并且在原字符串S中,L列是F列的前一个字符。我们根据该规则即可向前回溯得到原字符。回溯过程如下图所示:
最终我们可以得到字符串A2A1T0G0C0A0$,将该字符串反转就是我们的原字符S了。Python代码实现如下:
#!/usr/bin/env python3
# encoding: utf-8
"""
@author: burning
@software: pycharm
@file: decode_bwm.py
@time: 2020/5/25 20:54
@desc:
"""
def L_rank(rev_str):
total_num = {} # 字符出现的总次数,1开始
ranks = [] # 每个字符第几次出现,0开始
for i in rev_str: # 循环每个字符
if i not in total_num: # 当前字符未出现过
total_num[i] = 0
ranks.append(total_num[i]) # 添加当前字符第几次出现
total_num[i] += 1 # 当前字符出现次数加1
return total_num, ranks
def firstCol(total_num): # 获取F列
first = {}
total_c = 0
for c, count in sorted(total_num.items()): # 遍历排序后的L列
first[c] = (total_c, total_c+count) # 保存当前字符在F列中出现的位置范围
total_c += count # 下一个字符起始位置
return first
def reverse(str_L):
total_num, ranks = L_rank(str_L)
str_F = firstCol(total_num)
result="$" # 最终结果
start_row = 0 # 起始行
while str_L[start_row] != "$": # 非结束
c = str_L[start_row] # 当前字符
result = c + result # 回溯所得结果
start_row = str_F[c][0] + ranks[start_row]
return result
问答:
1、为什么通过转换后可以提高压缩水平?
答:因为转换后,大量相同字符将靠近,压缩时将会减少大量空间,例如:假设ATGC 在待压缩字符中大量存在,则通过该算法转换后,以TGC开头的字符串在排序后将会靠在一起,其最后的一个字符是A也会靠近在一起。
PS:还有部分是FM-index的原理,下次有机会再写了!
最新评论:
发表评论
电子邮件地址不会被公开。 必填项已用*标注