【Python实现】Smith-Waterman算法快速计算序列相似性

目录

• 1 SW算法介绍
• 2 ssw_aligner环境配置
• 3 项目安装地址
• 4 项目实例
• • 4.1Benchmark script:
  • 4.2 Input should be DNA type
• 5 参考文献

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)

5 参考文献

[1]生物信息学（2）——双序列比对之Smith-Waterman(SW)算法详解