比较序列相似性(sequence similarity)可以考虑用biopython或者emboss的几种比对方法。
1. Bio.pairwise2
主要用到SeqIO.parse读取,然后用Bio.pairwise2.align.globalxx比对并输出两个序列一样的比例。 如果用局部比对,可以用Bio.pairwise2.align.localxx.
from __future__ import division
from Bio import pairwise2 as pw2
from Bio import SeqIO
first_dict = SeqIO.to_dict(SeqIO.parse(open(first_fasta),'fasta')) # 直接转为字典格式
second_dict = SeqIO.to_dict(SeqIO.parse(open(second_fasta),'fasta'))
# 两个fasta文件中的序列两两比较:
for t in first_dict:
t_len = len(first_dict[t].seq)
for t2 in correspond[t]:
global_align = pw2.align.globalxx(first_dict[t].seq, second_dict[t2].seq)
matched = global_align[0][2]
percent_match = (matched / t_len) * 100
print(t + '\t' + t2 + '\t' + str(percent_match) + '\n')
2. Bio.Emboss.Applications
用了NeedleCommandline去比对,实测比上面的方法要快一点。不过都是python写的,又是基于DP,都不算很快。
from Bio.Emboss.Applications import NeedleCommandline
needle_cline = NeedleCommandline(asequence="test1.txt", bsequence="test2.txt",gapopen=10,gapextend=0.5,outfile='stdout')
out_data, err = needle_cline()
out_split = out_data.split("\n")
p = re.compile("\: (.*)\/")
print(int(p.search(out_split[24]).group(1).replace("%", "")))
3. needle
本质与上面的方法一样,不过这个是在shell中运行的。
# 安装
conda install -c bioconda emboss
# 运行
needle -asequence test1.txt -bsequence test2.txt -gapopen 10 -gapextend 0.5 -outfile aln.needle
grep 'Similarity' aln.needle | awk -F '[:| |/]' '{print $4}'