使用Biopython调用EMBOSS的needle

最近有一个需求,计算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,示例代码如下

1
2
3
4
5
6
7
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程序的路径。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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")

# build needle commnad
if needle_path is None:
needle_cmd = NeedleCommandline( asequence=aseq_file, bsequence=bseq_file,
gapopen=10, gapextend=0.5, outfile=needle_file)
else:
assert os.path.isfile(needle_path)
needle_cmd = NeedleCommandline( needle_path, asequence=aseq_file, bsequence=bseq_file,
gapopen=10, gapextend=0.5, outfile=needle_file)

# running needle_cmd
stdout, stderr = needle_cmd()
if verbose:
print(stdout + stderr, file = sys.stderr)

align = AlignIO.read(needle_file, "emboss")
# deleting the tempfile
os.unlink(aseq_file)
os.unlink(bseq_file)
os.unlink(needle_file)
return align

在写这个函数时,我在输出文件的文件名上折腾了很长时间,最初想使用原序列的ID,但ID包括一些奇奇怪怪的字符,例如(|:, 而NeedleCommandline不会自动在构建的命令中增加引号,使得运行时会出错。(你可以试试echo a(b会出现什么情况 )

后来,我使用了随机数作为文件名。这个方法在单个进程时问题不大,当时当我为了加速使用多进程时,却发现发现随机数重复了。进程A输出的文件A由于和进程B的输入文件名一样,导致被跑完的进程B删了。当然,解决方法也很简单,加上一些判断文件是否存在的ifelse语句即可。

不过,我还是选择调用tempfile标准库,使用mkstemp获取临时文件名,才解决了输出文件名的危机。

最终我就可以在循环中调用这个函数,然后将结果保存在列表中,用于后续处理。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
fasta = "/path/to/your/sequence.fasta"
seqs = [ fas for fas in SeqIO.parse(fasta, "fasta") ]
needle_path = "/opt/biosoft/EMBOSS-6.6.0/bin/needle"

align_res = []

for i in range(len(seqs)):
for j in range(i+1, len(seqs)):
align = run_needle(i, j, seqs , needle_path)
align_res.append(align )
# 以pickle保存python对象
outfile = open("pairwise_stats.txt", "w")
for align in align_res:
line = get_psa_stat(align)
line = "\t".join(line) + "\n"
outfile.write(line)

outfile.close()

根据我的计算,循环要运行3万多次,按照平均1s一个循环,只要500多分钟,小于9个小时,我只要在回去前运行这个程序(写完第一版大概是晚上10点),然后第二天早上过来收结果就行了。

但是,正如前文提及到的,「该程序会在多个进程运行时报错」,你也就知道我没有跑循环,而是选择了通过并行对循环进行加速。在一篇中,我将会讲一个,如何把一个简单的问题复杂化,最终又回到最初代码的故事