序列比对是生物信息学分析中的常见任务,包含局部比对和全局比对两大算法,局部比对最经典的代表是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,非常的方便。
|