In [3]:
from Bio import SeqIO, pairwise2
handle = open("../data/Reads.fa", "r")
records = list(SeqIO.parse(handle, "fasta"))
handle.close()
for record in records:
print record.id
print record.seq
In [4]:
seq1 = str(records[0].seq)
seq2 = str(records[1].seq)
alignments = pairwise2.align.globalxx(seq1, seq2)
print alignments
In [5]:
from Bio.pairwise2 import format_alignment
for aln in alignments : print(format_alignment(*aln))
这里出现了两个最佳比对，它们的Score都是44。 在每个最佳比对中，用竖线直观显示字符间对齐。 '-'表示空格， 例如，在第一个最佳比对中， read1的"AC-A"与read2的"ACCA"对齐，在第2个位置(位置从0开始)处read2比read1多了一个字符'C'。
我们管read2比read1在第2个位置处多了一个字符的情况，称为1个插入(insertion)，反之，称为删除(deletion)。 对齐的字符如果相同，称为match，不同称为mismatch。
insertion: '_' => 'A'
deletion： 'A' => '_'
match: 'A' => 'A'
mismatch: 'A' => 'C'
In [6]:
edit = [(i, j) for i, j in zip(alignments[0][0], alignments[0][1]) if i != j]
for e in edit: print e[0], '->', e[1]
In [42]:
S = 'ATCGAAGTG'
为了计算的方便，我们给s加一个'$'，用来标示字符串结尾。
In [43]:
S += '$'
SA(suffix array)是一个数组，它是以字典序的方式保存了S中所有后缀在S中的开始位置。
下面继续看例子，我们用Pos数组存放了S的所有后缀的开始位置。
In [44]:
Pos = range(len(S))
print Pos
以字典序的方式对Pos数组排序。
In [45]:
SA = sorted(Pos, key=lambda x:S[x:])
print SA
再看一下SA这些位置对应的后缀。
In [46]:
for i in SA: print S[SA[i]:]
放在一起对比一下。
In [47]:
print 'Index\tSA\tSuffix'
for i, p in enumerate(SA): print '{0}\t{1}\t{2}'.format(i, p, S[p:])
In [48]:
BWT = ''.join([S[x-1]for x in SA])
print BWT
再放一起看一下。
In [49]:
print 'Index\tSA\tBWT\tSuffix'
for i, p in enumerate(SA): print '{0}\t{1}\t{2}\t{3}'.format(i, p, S[p-1],S[p:])
发现了什么？
BWT[i]的字符是字符串s中SA[i]-1位置的字符，而Suffix[i]是S中SA[i]位置开始的后缀，也就是说其实BWT[i]是Suffix[i]的前一个字符。 为了更清楚的观察，我们这回用BWM矩阵，代表以字典序的方式排列S中的所有轮转(rotation)。
In [52]:
print 'Index\tSA\tBWT\tBWM'
for i, p in enumerate(SA): print '{0}\t{1}\t{2}\t{3}'.format(i, p, s[p-1],(s+s)[p:p+len(s)])
对比BWT和BWM发现，其实BWT就是BWM的最后一列。(想一想可以直接用SA来生成BMW)
In [58]:
firstCol = ''.join(sorted(BWT))
print firstCol
BWT作为BWM的最后一列，和firstCol有啥关系呢？
我们换一个方式来表示BWM，我们用字符在S上的位置表示每个rotation。
In [63]:
print 'Index\tSA\tBWT\tBWM\'\tBWM'
for i, p in enumerate(SA):
print '{0}\t{1}\t{2}\t{3}\t{4}'.format(i, p, s[p-1],''.join(map(str, range(len(S))[p:]+range(len(s))[:p])),(s+s)[p:p+len(s)])
对比BWM'与BWM发现以下的性质：
In [ ]: