BWT 算法和序列比对的基本实现

昨天晚上和今天抽空实现了Burrows Wheleer Tansform,并且尝试利用BWT，将短序列比对到长序列中。BWT的核心我觉得是要理解两个原则：

1. F序列的每个元素是下标对应的L元素的后一位。

2. 排序后，F中第一个A和L中第一个A是同一个A。（排序不改变相对位置），公共前缀不改变排序位置。

mapping 过程实现的非常基础，只能全序列不对，不能有gap。

#!/usr/bin/env python3
'''
burrows wheleer transform
核心是 排序后的F,L数组， L数组的某个元素是F数组对应的元素的前一位。
其次是，F数组的某个元素A在所有A中的相对位置不变。也就是F中的第一个A对应着L中的第一个A.
'''import os
import sys
from collections import Counterdef BWT(s, end):'''输出1：L string输出2：L string 每个字符的原来的位置'''# 记录原序列中字符出现位置C = {i: [] for i in set(s)}  # 初始化L列的元素相对位置for i in range(len(s)):  # 分别给每个元素标记出现的位置C[s[i]].append(i)s = s + end  # 添加最后的符号M = []  # 初始化输出count = len(s)  # 计算字符串长度# 设置偏移矩阵while True:M.append(s)s = s[len(s)-1]+s[0:len(s)-1]count -= 1if count <= 0:breakN = [[v, k] for k, v in enumerate(M[::-1])] # 标记位置信息N_sort = sorted(N, key=lambda x: x[0])  #字典排序C = {i: [] for i in set(s)}  # 初始化L列的元素相对位置L = []for k, v in N_sort: # 以A为例，Lstring中 从头开始每一个A在原来字符串的位置 L.append(k[-1])C[k[-1]].append(v)return(["".join(L), C])def reverseBWT(s, end):''''''B = Counter(s)  # 压缩的F列C = {i: [] for i in set(B)}  # 初始化L列的元素相对位置for i in range(len(s)):  # 分别给每个元素标记出现的位置C[s[i]].append(i)out = ""  # 初始化输出arrow = 1while True:Lbase = s[arrow-1]  # 取出对应的字符out = out + Lbase  # 添加到输出if Lbase == end:  # 如果遇到终止符，则退出循环，并输出break# 从数组C中查找对应字符出现的位置，并且枚举，这样可以得到字符所在位置对应的相同字符的偏移量。for i, j in enumerate(C[Lbase]):if j == arrow-1:pianyi = i+1  # 这里应该是对应的L中的序号breakarrow = CheckF(B, Lbase, pianyi)return(out[::-1])def simple_mapping(ref, seed):s, rank = BWT(ref, "#")B = Counter(s)  # 压缩的F列C = {i: [] for i in set(s)}  # 初始化L列的元素相对位置for i in range(len(s)):  # 分别给每个元素标记出现的位置C[s[i]].append(i)seed = seed[::-1]print(ref)for mark in range(len(C[seed[0]])):arrow = C[seed[0]][mark]  # mark:0,1,2 编号，0开始# arrow: L 中的 12,34,65... 实际位置out = ""count = 0while True:Lbase = s[arrow]if Lbase != seed[count] or Lbase == '#':breakout = out + Lbasefor i, v in enumerate(C[Lbase]):if v == arrow:pianyi = i + 1breakarrow = CheckF_0base(B, Lbase, pianyi)count += 1if count > len(seed)-1:left = rank[seed[0]][mark]out = "*"*(left-len(out)+1) + out[::-1] + "*"*(len(s)-left-2)print(out)breakdef CheckF(mtx, base, order):out_order = 0for k in sorted(mtx.keys()):  # 从F数组中根据字符，以及偏移量，计算出下一个坐标 P = 字符前所有的字符的数目+偏移量if base != k:out_order += mtx[k]else:out_order += orderbreakreturn(out_order)def CheckF_0base(mtx, base, order):out_order = 0for k in sorted(mtx.keys()):  # 从F数组中根据字符，以及偏移量，计算出下一个坐标 P = 字符前所有的字符的数目+偏移量if base != k:out_order += mtx[k]else:out_order += orderbreakreturn(out_order-1)if __name__ == '__main__':x = "actgagctttagcgtagctttaggagagcttcctagctacgtatcgagcgggcatctatc"seed = "ga"print("origin string is :" + x)print("L string      is :" + BWT(x, '#')[0])print("seed string   is :" + seed)print("mapped seq    is :")simple_mapping(x, seed)

结果：

BWT 算法和序列比对的基本实现相关推荐

序列比对（27）BWT算法
本文介绍了BWT算法. bwa是目前最流行的二代测序比对工具,其中就用到了BWT算法.BWT(Burrows-Wheeler Transform)算法是一种数据转换算法,它将一个字符串中的相似字符放在 ...
BWT算法解析及Java语言实现
BWT算法将原来的文本转换为一个相似的文本,转换后使得相同的字符位置连续或者相邻,之后可以使用其他技术如:Move-to-fronttransform 和游程编码进行文本压缩. BWT原理: 1. ...
北京大学生物信息学-第五周-新一代测序(NGS) 回帖 BWT算法
新一代测序 Read: A short DNA fragment which is read out by sequencer. 读:由测序仪读出的短DNA片段. DNA序列+质量信息->FAS ...
Algs4-1.5.1使用quick-find算法处理序列
1.5.1使用quick-find算法处理序列9-0 3-4 5-8 7-2 2-1 5-7 0-3 4-2.对于输入的每一对整数,给出id[]数组的内容和访问数组的次数. 答: public cla ...
Python使用爬山算法寻找序列“最大值”
爬山算法是人工智能算法的一种,特点在于局部择优,所以不一定能够得到全局最优解,尽管效率比较高.使用爬山算法寻找序列最大值的思路是:在能看得到的局部范围内寻找最大值,如果当前元素已经是最大值就结束,如果 ...
手把手教你深度学习强大算法进行序列学习(附Python代码)
作者:NSS 翻译:陈之炎校对:丁楠雅本文共3200字,建议阅读10分钟. 本文将教你使用做紧致预测树的算法来进行序列学习. 概述序列学习是近年来深度学习的热点之一.从推荐系统到语音识别再到自然 ...
hdu 4358(莫队算法+dfs序列)
题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=4358 解题思路:用dfs求出整棵树的dfs序列,这样以u为根节点的子树就转化到相对应的区间上了.由于是 ...
viterbi算法_序列比对(十四)——viterbi算法和后验解码的比较
原创: hxj7 本文比较了viterbi算法求解最可能路径以及后验解码这两种不同的解码方法. 前文<序列比对(十)viterbi算法求解最可能路径>介绍了用viterbi算法求解最可能路 ...
蒙特卡洛粒子滤波定位算法_序列蒙特卡洛(SMC)与粒子滤波
上一节我们讲了蒙特卡洛方法,如果对蒙特卡洛方法不了解的朋友们可以看我之前写的这一篇文章. Shawn:一文看懂蒙特卡洛采样方法zhuanlan.zhihu.com 今天这篇文章我们来讲一讲针对时间序 ...

BWT 算法和序列比对的基本实现

BWT 算法和序列比对的基本实现相关推荐

最新文章

热门文章