最近有一个需求，计算250条序列两两之间的相似度，之后根据相似度继续进行序列过滤，即A和B如果相似度高于一个阈值，那么就把B去掉。
我最初想到的就是在Python里面调用EMBOSS的needle，准备使用Python的os或者subprocess模块，但是突然想到之前看WGDI源代码时，发现作者时通过biopython的API进行MAFFT的调用，于是就决定改用Biopython.
通过查找文档，我找到了Biopython对这部分内容的介绍，见 http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec97,示例代码如下
from Bio.Emboss.Applications import NeedleCommandline
needle_cline = NeedleCommandline(asequence="alpha.faa",
bsequence="beta.faa",
gapopen=10, gapextend=0.5,
outfile="needle.txt")
stdout, stderr = needle_cline()
align = AlignIO.read("needle.txt", "emboss")
但是这个代码距离我最终代码还有很长的路要走，因为我需要先将序列写出到文件，然后调用命令行，然后再读取输出文件。我的问题是，能不能省掉这几次I/O操作。在Emboss文档的下面，我发现Biopython提供了两个模块，分别是Bio.pairwise2和Bio.Align.PairwiseAligner, 可以直接使用Biopython的.SeqRecord对象。按照文档描述，pairwise2和PairwiseAligner能够达到和调用EMBOSS相同的速度，经过我测试，发现两者速度的确差别不大，这说明直接在Python里面进行配对序列联配并不比先输出运行所需文件，然后读取运行生成文件这个流程快，所以经过一番折腾之后，我还是决定使用NeedleCommandline。
为了方便后续调用，我写了一个函数，函数的前3个参数是用来从SeqRecord列表中提取序列用于写出，第4个参数用来提供needle程序的路径。
def run_needle(aseq_pos, bseq_pos, seq, needle_path, verbose=False):
"""run EMBOSS needle retrun align
"""
# get the Bio object
aseq = seq[aseq_pos]
bseq = seq[bseq_pos]
# create tempfile for write fasta and needle result
aseq_file = mkstemp(suffix=".fas")[1]
bseq_file = mkstemp(suffix=".fas")[1]
needle_file = mkstemp(suffix=".needle")[1]
# write the sequence
SeqIO.write(aseq, aseq_file, "fasta")
SeqIO.write(bseq, bseq_file, "fasta"
