Biopython 序列比对

序列比对 是将两个或多个序列（DNA、RNA或蛋白质序列）按照特定的顺序排列，以确定它们之间的相似区域的过程。

识别相似区域使我们能够推断出很多信息，如物种之间有哪些性状是保守的，不同物种在遗传上有多接近，物种是如何进化的，等等。Biopython为序列比对提供了广泛的支持。

让我们在本章中学习Biopython提供的一些重要功能 –

解析序列比对

Biopython提供了一个模块，Bio.AlignIO来读写序列排列。在生物信息学中，有很多格式可以用来指定序列排列数据，类似于早期学习的序列数据。Bio.AlignIO提供的API与Bio.SeqIO类似，只是Bio.SeqIO是在序列数据上工作，而Bio.AlignIO是在序列排列数据上工作。

在开始学习之前，让我们从网上下载一个序列比对样本文件。

要下载该样本文件，请按照以下步骤进行

第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 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
>>>

我们还可以检查排列中可用的序列（SeqRecord），以及下面的-

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

多重排列

一般来说，大多数的序列比对文件都包含单一的比对数据，使用读的方法来解析它就足够了。在多序列对齐的概念中，两个或更多的序列被比较，以获得它们之间的最佳子序列匹配，并在一个文件中产生多序列对齐。

如果输入的序列排列格式包含一个以上的序列排列，那么我们需要使用解析方法，而不是如下所述的读取方法

>>> from Bio import AlignIO 
>>> alignments = AlignIO.parse(open("PF18225_seed.txt"), "stockholm") 
>>> print(alignments) 
<generator object parse at 0x000001CD1C7E0360> 
>>> 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步

与seq1和seq2一起调用pairwise2.align.globalxx方法，用下面这行代码找到排列组合

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

在这里， globalxx 方法执行实际工作，在给定的序列中找到所有可能的最佳排列。实际上，Bio.pairwise2提供了相当多的方法，这些方法遵循以下惯例来寻找不同情况下的排列组合。

<sequence alignment type>XY

这里，序列排列类型指的是排列类型，可以是 全局的 或 局部的。全局 类型是通过考虑整个序列来寻找序列排列。局部类型是通过查看给定序列的子集来寻找序列排列。这将是繁琐的，但能更好地了解给定序列之间的相似性。

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步

在可迭代的alignments对象上循环，得到每个单独的对齐对象并打印出来。

>>> 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 needle 和 water

让我们在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，输入文件opuntia.fasta可在Biopython软件包中找到。 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 columns
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273285|gb|AF191659.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273284|gb|AF191658.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273287|gb|AF191661.1|AF191
TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273286|gb|AF191660.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273290|gb|AF191664.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273289|gb|AF191663.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA
gi|6273291|gb|AF191665.1|AF191
>>>