1 SW算法介绍
Smith-Waterman 算法是由 Temple F. Smith 和 Michael S. Waterman 两人在 1981 年提出来的,是 Needleman-Wunsch 算法的改良版,通过算法的比对,能获取到局部最优解。具体的算法这里不做过多解释,详见参考文献[1]。本文关注的主要问题是,计算速度的问题。尽管网上有很多SW算法的实现,有Java、Python和R包。然而,当在计算大量的序列相似性时,以上代码的实现会非常的慢,甚至需要月余。如果使用C语言实现代码,效果非常明显但是可扩展性不高。考虑到以上两个因素,Python下的Cython库是一个非常好的实现方式,可以实现Python和C语言的混合编程。
2 ssw_aligner环境配置
### Dependencies
- [numpy==1.12.0](http://www.numpy.org/)
- [Cython==0.28.3](https://cython.org/)
3 项目安装地址
https://github.com/mbreese/swalign
4 项目实例
from ssw_aligner import local_pairwise_align_ssw
query_seq = 'TTTTTAAAAA'
target_seq = 'GGGGTTTT'
alignment = local_pairwise_align_ssw(query_seq,
target_seq,
gap_open_penalty=11,
gap_extend_penalty=1,
match_score=2,
mismatch_score=-3)
# get score
alignment.optimal_alignment_score
# get query start, end
alignment.query_begin
alignment.query_end
# get target start, end
alignment.target_begin
alignment.target_end_optimal
# get aligned sequence
alignment.aligned_query_sequence
alignment.aligned_target_sequence
# get cigar infomation
alignment.cigar
# check whether the index starts from 0 or not
alignment.is_zero_based()
# make the index start from n(0 or 1)
alignment.set_zero_based(0) # start from 0
alignment.set_zero_based(1) # start from 1
4.1Benchmark script:
import random
import time
from skbio import DNA
import skbio
import swalign
import ssw_aligner
match = 2
mismatch = -1
scoring = swalign.NucleotideScoringMatrix(match, mismatch)
sw = swalign.LocalAlignment(scoring)
bases = ['A', 'T', 'C', 'G']
def generate_gene(length):
return ''.join([random.choice(bases) for i in range(0, length)])
def benchmark(align_func):
start = time.time()
for i in range(0, 100):
for seq_length in range(100, 2000, 500):
seq1, seq2 = generate_gene(seq_length), generate_gene(seq_length)
align_func(seq1, seq2)
return time.time() - start
4.2 Input should be DNA type
def benchmark_skbio(align_func):
start = time.time()
for i in range(0, 100):
for seq_length in range(100, 2000, 500):
seq1, seq2 = generate_gene(seq_length), generate_gene(seq_length)
align_func(DNA(seq1), DNA(seq2))
return time.time() - start
print('ssw_aligner')
ssw_aligner_time = benchmark(ssw_aligner.local_pairwise_align_ssw)
print(ssw_aligner_time)
print('skbio')
skbio_time = benchmark_skbio(skbio.alignment.local_pairwise_align_ssw)
print(skbio_time)
print('swalign')
swalign_time = benchmark(sw.align)
print(swalign_time)