欢迎关注”生信修炼手册”!
序列比对是生物信息学分析中的常见任务，包含局部比对和全局比对两大算法，局部比对最经典的代表是blast, 全局比对则用于多序列比对。在biopython中，支持对序列比对的结果进行读写，解析，以及运行序列比对的程序。
首先来看下多序列比对，多序列比对的软件较多，比如clustalw, muscle, mafft等，输出结果的格式也很多，比如clustal, fasta, phylip等。在biopython中，为不同格式，不同软件提供了统一的接口，方便我们的使用
1. 读取多序列比对结果
通过Bio.AlignIO模块来对多序列比对结果进行读写，其中的parse方法用于从文件句柄中读取多序列比对的内容，用法如下
-
>>> from Bio
import AlignIO
-
>>> alignment = AlignIO.parse(
'clustal.out',
'clustal')
-
>>>
print(alignment)
-
<generator object parse at
0x0928C300>
-
>>>
for i in alignment:
-
...
print(i.id)
-
...
该方法的返回值是一个迭代器，每次迭代，返回的是一个SeqRecord对象。
2. 输出多序列比对结果
通过write方法将多序列比对的结果输出到文件中，可以指定输出文件的格式，用法如下
-
>>> alignments = AlignIO.parse(
"aln.fasta",
"fasta")
-
>>> AlignIO.write(alignments,
"aln.clustal",
"clustal")
和Bio.SeqIO相同，针对格式转换，也体用了convert方法，用法如下
>>> count = AlignIO.convert("aln.fasta", "fasta", "align.clustal", "clustal")3. 运行多序列比对程序
为了简化调用，在Bio.Applicaitons模块中，提供了各种应有的调用接口。对于多序列比对结果，默认调用位于本地PATH环境变量下的可执行程序，来执行对应的命令，以clustalw为例，用法如下
-
>>> from Bio.Align.Applications
import ClustalwCommandline
-
>>> cline = ClustalwCommandline(
"clustalw2", infile=
"input.fasta")
第一个参数指定可执行程序，如果可执行程序位于PATH变量下，指定可执行程序的名称即可，否则需要指定可执行程序的完整路径。clustalw会根据输入文件的名称，自动确定输出文件的名字。当然，也可以通过参数指定输出文件的名字。
Bio.Applicaitons模块通过subprocess来调用程序，我们可以借此来读取程序的标准输出和标准错误流信息。比如以muscle为例，通过直接读取标准输出，避免了创建临时文件，示例如下
-
>>>
import subprocess
-
>>> from Bio.Align.Applications
import MuscleCommandline
-
>>> from Bio
import AlignIO
-
>>> muscle_cline = MuscleCommandline(input=
"opuntia.fasta")
-
>>> child = subprocess.Popen(str(muscle_cline), stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True)
-
>>> align = AlignIO.read(child.stdout,
"fasta")
对于局部比对而言，biopython可以运行blast并解析其输出
1. 运行blast
支持联网运行和本地运行两种模式，联网运行时调用NCBI网站的blast程序，用法如下
-
# 传统的文件读取, 适合fasta格式
-
>>> from Bio.Blast
import NCBIWWW
-
>>> fasta_string = open(
"input.fasta").read()
-
>>> result_handle = NCBIWWW.qblast(
"blastn",
"nt", fasta_string)
-
# Bio.SeqIO读取，适合fasta,genebank等格式
-
>>> record = SeqIO.read(
"input.fasta", format=
"fasta")
-
>>> result_handle = NCBIWWW.qblast(
"blastn",
"nt", record.format(
'fasta'))
在线运行只需要我们提供查询序列即可，用的数据库是NCBI的公共数据库，而本地运行则要求我们在本地安装好blast可执行程序，构建本地版本的blast数据库才行，本地运行的用法如下
-
>>> from Bio.Blast.Applications
import NcbiblastxCommandline
-
>>> blastx_cline = NcbiblastxCommandline(query=
"query.fasta", db=
"nr", evalue=
0.001, outfmt=
5, out=
"output.xml")
-
>>> stdout, stderr = blastx_cline()
2. 解析blast的输出
biopython中blast默认的输出格式为xml, 解析其输出的用法如下
-
>>> from Bio.Blast
import NCBIXML
-
>>> blast_records = NCBIXML.parse(result_handle)
-
>>> E_VALUE_THRESH =
0.001
-
>>>
for blast_record in blast_records:
-
...
for alignment in blast_record.alignments:
-
...
for hsp in alignment.hsps:
-
...
if hsp.expect < E_VALUE_THRESH:
-
...
print
'****Alignment****'
-
...
print
'sequence:', alignment.title
-
...
print
'length:', alignment.length
-
...
print
'e value:', hsp.expect
-
...
print hsp.query[
0:
75] +
'...'
-
...
print hsp.match[
0:
75] +
'...'
-
...
print hsp.sbjct[
0:
75] +
'...'
对于序列比对结果的运行和解析，通过biopython可以很好的将其整合到python生态中，对于用python构建一套完整的pipeline,非常的方便。
·end·
—如果喜欢，快分享给你的朋友们吧—
原创不易，欢迎收藏，点赞，转发！生信知识浩瀚如海，在生信学习的道路上，让我们一起并肩作战！
本公众号深耕耘生信领域多年，具有丰富的数据分析经验，致力于提供真正有价值的数据分析服务，擅长个性化分析，欢迎有需要的老师和同学前来咨询。
更多精彩
写在最后
转发本文至朋友圈，后台私信截图即可加入生信交流群，和小伙伴一起学习交流。
扫描下方二维码，关注我们，解锁更多精彩内容！
一个只分享干货的
生信公众号
转载：https://blog.csdn.net/weixin_43569478/article/details/112001319