Biopython - 序列比对

序列比对是按特定顺序排列两个或多个序列(DNA，RNA或蛋白质序列)的过程，以确定它们之间的相似区域.

识别相似区域使我们能够推断出大量信息，例如物种之间保守的性状，遗传上不同物种的接近程度，物种如何进化等.Biopython为序列比对提供了广泛的支持.

让我们在本章中了解Biopython提供的一些重要特性 :

解析序列比对

Biopython提供了模块，Bio.AlignIO读取和写入序列比对.在生物信息学中，有许多格式可用于指定与早期学习的序列数据类似的序列比对数据. Bio.AlignIO提供类似于Bio.SeqIO的API，除了Bio.SeqIO处理序列数据，Bio.AlignIO处理序列比对数据.

在开始学习之前，让我们从Internet下载样本序列比对文件.

要下载样本文件，请按照以下步骤 : 去;

步骤1 : 打开您喜欢的浏览器并转到 http://pfam.xfam.org/family/browse 网站.它将按字母顺序显示所有Pfam系列.

第2步 : 选择具有较少种子值的任何一个家庭.它包含最少的数据，使我们能够轻松地进行对齐.在这里，我们选择/点击了PF18225并打开了 http://pfam.xfam.org/family/PF18225 并显示有关它的完整详细信息，包括序列比对.

步骤3 : 转到对齐部分并以斯德哥尔摩格式(PF18225_seed.txt)下载序列比对文件.

让我们尝试使用Bio.AlignIO读取下载的序列比对文件，如下所示 :

导入Bio.AlignIO模块

>>> from Bio import AlignIO

使用read方法读取对齐. read方法用于读取给定文件中可用的单个对齐数据.如果给定文件包含许多对齐，我们可以使用parse方法. parse方法返回可迭代的对齐对象，类似于Bio.SeqIO模块中的parse方法.

>>> alignment = AlignIO.read(open("PF18225_seed.txt"), "stockholm")

打印对齐对象.

>>> print(alignment)SingleLetterAlphabet() alignment with 6 rows and 65 columnsMQNTPAERLPAIIEKAKSKHDINVWLLDRQGRDLLEQRVPAKVA...EGP B7RZ31_9GAMM/59-123 AKQRGIAGLEEWLHRLDHSEAIPIFLIDEAGKDLLEREVPADIT...KKP A0A0C3NPG9_9PROT/58-119 ARRHGQEYFQQWLERQPKKVKEQVFAVDQFGRELLGRPLPEDMA...KKP A0A143HL37_9GAMM/57-121 TRRHGPESFRFWLERQPVEARDRIYAIDRSGAEILDRPIPRGMA...NKP A0A0X3UC67_9GAMM/57-121 AINRNTQQLTQDLRAMPNWSLRFVYIVDRNNQDLLKRPLPPGIM...NRK B3PFT7_CELJU/62-126 AVNATEREFTERIRTLPHWARRNVFVLDSQGFEIFDRELPSPVA...NRT K4KEM7_SIMAS/61-125>>>

我们还可以检查对齐中可用的序列(SeqRecord)以及低于和低于;

>>> for align in alignment: ... print(align.seq) ... MQNTPAERLPAIIEKAKSKHDINVWLLDRQGRDLLEQRVPAKVATVANQLRGRKRRAFARHREGP AKQRGIAGLEEWLHRLDHSEAIPIFLIDEAGKDLLEREVPADITA---RLDRRREHGEHGVRKKP ARRHGQEYFQQWLERQPKKVKEQVFAVDQFGRELLGRPLPEDMAPMLIALNYRNRESHAQVDKKP TRRHGPESFRFWLERQPVEARDRIYAIDRSGAEILDRPIPRGMAPLFKVLSFRNREDQGLVNNKP AINRNTQQLTQDLRAMPNWSLRFVYIVDRNNQDLLKRPLPPGIMVLAPRLTAKHPYDKVQDRNRK AVNATEREFTERIRTLPHWARRNVFVLDSQGFEIFDRELPSPVADLMRKLDLDRPFKKLERKNRT >>>

多个对齐

通常，大多数序列比对文件都包含单个比对数据，它足以使用读取方法来解析它.在多序列比对概念中，比较两个或更多序列以获得它们之间的最佳子序列匹配，并在单个文件中导致多序列比对.

如果输入序列比对格式包含多个序列比对，然后我们需要使用 parse 方法而不是 read 方法，如下所示 :

>>> from Bio import AlignIO >>> alignments = AlignIO.parse(open("PF18225_seed.txt"), "stockholm") >>> print(alignments)  >>> for alignment in alignments: ... print(alignment) ... SingleLetterAlphabet() alignment with 6 rows and 65 columns MQNTPAERLPAIIEKAKSKHDINVWLLDRQGRDLLEQRVPAKVA...EGP B7RZ31_9GAMM/59-123 AKQRGIAGLEEWLHRLDHSEAIPIFLIDEAGKDLLEREVPADIT...KKP A0A0C3NPG9_9PROT/58-119 ARRHGQEYFQQWLERQPKKVKEQVFAVDQFGRELLGRPLPEDMA...KKP A0A143HL37_9GAMM/57-121 TRRHGPESFRFWLERQPVEARDRIYAIDRSGAEILDRPIPRGMA...NKP A0A0X3UC67_9GAMM/57-121 AINRNTQQLTQDLRAMPNWSLRFVYIVDRNNQDLLKRPLPPGIM...NRK B3PFT7_CELJU/62-126 AVNATEREFTERIRTLPHWARRNVFVLDSQGFEIFDRELPSPVA...NRT K4KEM7_SIMAS/61-125>>>

这里，parse方法返回可迭代的对齐对象，可以迭代它以获得实际的对齐.

成对序列对齐

成对序列比对一次只比较两个序列，并提供最佳的序列比对. 成对易于理解，并且可以从结果序列比对中推断出来.

Biopython提供了一个特殊的模块， Bio.pairwise2 使用成对方法识别比对序列. Biopython应用最佳算法来找到比对序列，它与其他软件相同.

让我们编写一个例子，使用成对模块找到两个简单和假设序列的序列比对.这将有助于我们理解序列比对的概念以及如何使用Biopython对其进行编程.

步骤1

导入模块 pairwise2 使用下面给出的命令 :

>>> from Bio import pairwise2

步骤2

创建两个序列，seq1和seq2 :

>>> from Bio.Seq import Seq >>> seq1 = Seq("ACCGGT") >>> seq2 = Seq("ACGT")

步骤3

调用方法pairwise2.align.globalxx以及seq1和seq2使用下面的代码行找到对齐 :

 >>> alignments = pairwise2.align.globalxx(seq1，seq2)

此处， globalxx 方法执行实际工作并找到最佳效果给定序列中可能的比对.实际上，Bio.pairwise2提供了一系列遵循以下惯例的方法，以在不同情况下找到对齐.

XY

这里，序列比对类型是指对齐类型，可以是 global 或 local.全局类型通过考虑整个序列来发现序列比对.本地类型也通过查看给定序列的子集来寻找序列比对.这将是乏味的，但更好地了解给定序列之间的相似性.

• X指的是匹配分数.可能的值是x(完全匹配)，m(基于相同字符的分数)，d(用户提供的带字符和匹配分数的字典)，最后是c(用户定义的函数，用于提供自定义评分算法).

• Y指间隙罚款.可能的值是x(无间隙罚分)，s(两个序列的相同惩罚)，d(每个序列的不同惩罚)，最后c(用户定义的函数提供自定义间隙罚分)

因此，localds也是一个有效的方法，它使用局部对齐技术找到序列比对，用户提供的匹配字典和用户提供的两个序列的空位罚分.

 >>> test_alignments = pairwise2.align.localds(seq1，seq2，blosum62，-10，-1)

这里，blosum62指的是pairwise2模块中可用的字典提供比赛得分. -10表示缺口开放罚分，-1表示缺口延长罚金.

步骤4

循环遍历可迭代路线对象并得到每个人对齐对象并打印它.

>>> for alignment in alignments: ... print(alignment) ... ('ACCGGT', 'A-C-GT', 4.0, 0, 6) ('ACCGGT', 'AC--GT', 4.0, 0, 6) ('ACCGGT', 'A-CG-T', 4.0, 0, 6) ('ACCGGT', 'AC-G-T', 4.0, 0, 6)

步骤5

Bio.pairwise2模块提供格式方法，format_alignment更好地可视化结果 :

>>> from Bio.pairwise2 import format_alignment >>> alignments = pairwise2.align.globalxx(seq1, seq2) >>> for alignment in alignments: ... print(format_alignment(*alignment)) ...ACCGGT | | || A-C-GT    Score=4    ACCGGT || || AC--GT    Score=4 ACCGGT | || | A-CG-T    Score=4 ACCGGT || | | AC-G-T    Score=4>>>

Biopython还提供另一个模块进行序列比对，Align.这个模块提供了一组不同的API来简单地设置参数，如算法，模式，匹配分数，空位罚分等.简单查看Align对象如下&&;

>>> from Bio import Align>>> aligner = Align.PairwiseAligner()>>> print(aligner)Pairwise sequence aligner with parameters   match score: 1.000000   mismatch score: 0.000000   target open gap score: 0.000000   target extend gap score: 0.000000   target left open gap score: 0.000000   target left extend gap score: 0.000000   target right open gap score: 0.000000   target right extend gap score: 0.000000   query open gap score: 0.000000   query extend gap score: 0.000000   query left open gap score: 0.000000   query left extend gap score: 0.000000   query right open gap score: 0.000000   query right extend gap score: 0.000000   mode: global>>>

支持序列比对工具

Biopython通过Bio.Align.Applications提供了许多序列比对工具的接口模块.一些工具列在下面和下面;

• ClustalW

• MUSCLE

• EMBOSS针和水

让我们在Biopython中编写一个简单的例子，通过最流行的对齐工具ClustalW创建序列对齐.

第1步 : 从 http://www.clustal.org/download/current/下载Clustalw程序并安装它.另外，使用"clustal"安装路径更新系统PATH.

步骤2 : 从模块Bio.Align.Applications导入ClustalwCommanLine.

>>> from Bio.Align.Applications import ClustalwCommandline

步骤3 : 通过使用输入文件调用ClustalwCommanLine来设置cmd，在Biopython包中提供opuntia.fasta.  https://raw.githubusercontent.com/biopython/biopython/master/Doc/examples/opuntia.fasta

 >>> cmd = ClustalwCommandline("clustalw2"， infile ="/path/to/biopython/sample/opuntia.fasta")>>> print(cmd) clustalw2 -infile = fasta/opuntia.fasta

第4步 : 调用cmd()将运行clustalw命令并给出生成的
对齐文件的输出，opuntia.aln.

 > >> stdout，stderr = cmd()

第5步 : 阅读并打印对齐文件，如下所示;

>>> from Bio import AlignIO>>> align = AlignIO.read("/path/to/biopython/sample/opuntia.aln", "clustal")>>> print(align)SingleLetterAlphabet() alignment with 7 rows and 906 columnsTATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGAgi|6273285|gb|AF191659.1|AF191TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGAgi|6273284|gb|AF191658.1|AF191TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGAgi|6273287|gb|AF191661.1|AF191TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGAgi|6273286|gb|AF191660.1|AF191TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGAgi|6273290|gb|AF191664.1|AF191TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGAgi|6273289|gb|AF191663.1|AF191TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGAgi|6273291|gb|AF191665.1|AF191>>>