BioPython教程翻译：6.1分析或读取序列比对，6.1  Parsing or Reading Sequence Alignments

6.1  分析或读取序列比对

我们有两个用于读取序列比对的函数 ， Bio.AlignIO.read()和Bio.AlignIO.parse() ，它们与Bio.SeqIO 中的对应函数的分工是一样的 ，一个用来读取只有一条比对的文件，另一个用来读取有多条比对的文件。

使用Bio.AlignIO.parse()的话，会返回一个 迭代器 ，它指向 MultipleSeqAlignment 对象 。迭代器通常是在一个for 循环中使用的 。比如以下例子：你有很多不同的比对 ，它们包括来自PHYLIP 工具 seqboot 的重采样比对 、或者来自EMBOSS 工具 water 或 needle 或Bill Pearson 的FASTA 工具的多个成对的比对 。

然而，很多时候 ，你所研究的文件中只有一条比对。在这种情况下 ，你应当使用 Bio.AlignIO.read() 函数 ，它返回单个 MultipleSeqAlignment 对象 。

这2个函数都有2个必选参数：

1. 1.第一个参数是一个用来读取数据的把柄 （ handle ），典型地 ，是一个打开的文件（参见 20.1 小节），或者是一个文件名。

2. 2.第二个参数是一个小写的字符串，指定比对的格式。就像在 Bio.SeqIO 中一样 ，我们不会尝试为你猜测那个格式！参见 http://biopython.org/wiki/AlignIO ，那里有 被 支持的格式的完整列表 。

另外还有一个可选的seq_count 参数，我们将在下面的 6.1.3 小节中说明它，它是用来对付那些有歧义的可能包含多条比对信息的文件格式的 。

Biopython 1.49中还引入咯另一个可选参数alphabet ，它允许你指定一个预期的字符集 （alphabet）。这个选项会有用的 ，因为很多比对文件格式中 都没有显式指明它们的序列是RNA、DNA 还是蛋白质序列–那就意味着Bio.AlignIO 会使用默认的通用字符集。

6.1.1  单条比对信息

举个例子，看看下面带有丰富注释（annotation rich）的蛋白质比对信息，以PFAM 或Stockholm 文件格式存储的：

# STOCKHOLM 1.0

#=GS COATB_BPIKE/30-81  AC P03620.1

#=GS COATB_BPIKE/30-81  DR PDB; 1ifl ; 1-52;

#=GS Q9T0Q8_BPIKE/1-52  AC Q9T0Q8.1

#=GS COATB_BPI22/32-83  AC P15416.1

#=GS COATB_BPM13/24-72  AC P69541.1

#=GS COATB_BPM13/24-72  DR PDB; 2cpb ; 1-49;

#=GS COATB_BPM13/24-72  DR PDB; 2cps ; 1-49;

#=GS COATB_BPZJ2/1-49   AC P03618.1

#=GS Q9T0Q9_BPFD/1-49   AC Q9T0Q9.1

#=GS Q9T0Q9_BPFD/1-49   DR PDB; 1nh4 A; 1-49;

#=GS COATB_BPIF1/22-73  AC P03619.2

#=GS COATB_BPIF1/22-73  DR PDB; 1ifk ; 1-50;

COATB_BPIKE/30-81             AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA

#=GR COATB_BPIKE/30-81  SS    -HHHHHHHHHHHHHH--HHHHHHHH--HHHHHHHHHHHHHHHHHHHHH----

Q9T0Q8_BPIKE/1-52             AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA

COATB_BPI22/32-83             DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA

COATB_BPM13/24-72             AEGDDP...AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA

#=GR COATB_BPM13/24-72  SS    ---S-T...CHCHHHHCCCCTCCCTTCHHHHHHHHHHHHHHHHHHHHCTT--

COATB_BPZJ2/1-49              AEGDDP...AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA

Q9T0Q9_BPFD/1-49              AEGDDP...AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA

#=GR Q9T0Q9_BPFD/1-49   SS    ------...-HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH--

COATB_BPIF1/22-73             FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA

#=GR COATB_BPIF1/22-73  SS    XX-HHHH--HHHHHH--HHHHHHH--HHHHHHHHHHHHHHHHHHHHHHH---

#=GC SS_cons                  XHHHHHHHHHHHHHHHCHHHHHHHHCHHHHHHHHHHHHHHHHHHHHHHHC--

#=GC seq_cons                 AEssss...AptAhDSLpspAT-hIu.sWshVsslVsAsluIKLFKKFsSKA

//

这是Phage_Coat_Gp8 (PF05371)那个PFAM 条目的种子（seed）比对，是从现已过时的PFAM 版本 http://pfam.sanger.ac.uk/ 上下载回来的 。我们可以用下面的代码来载入这个文件 （假设它已经保存到当前工作目录中 ，文件名为 “PF05371_seed.sth”）：

>>> from Bio import AlignIO

>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")

以下代码会输入这个比对信息的概要内容：

>>> print alignment

SingleLetterAlphabet() alignment with 7 rows and 52 columns

AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRL...SKA COATB_BPIKE/30-81

AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKL...SRA Q9T0Q8_BPIKE/1-52

DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRL...SKA COATB_BPI22/32-83

AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPM13/24-72

AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA COATB_BPZJ2/1-49

AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKL...SKA Q9T0Q9_BPFD/1-49

FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKL...SRA COATB_BPIF1/22-73

你应该会注意到，在上面的输出内容中，那些序列都 被省略咯一部分。我们可以自己写代码来输出这些东西，把那些行当成 SeqRecord 对象来遍历 ：

>>> from Bio import AlignIO

>>> alignment = AlignIO.read("PF05371_seed.sth", "stockholm")

>>> print "Alignment length %i" % alignment.get_alignment_length()

Alignment length 52

>>> for record in alignment:

...     print "%s - %s" % (record.seq, record.id)

AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA - COATB_BPIKE/30-81

AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA - Q9T0Q8_BPIKE/1-52

DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA - COATB_BPI22/32-83

AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - COATB_BPM13/24-72

AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA - COATB_BPZJ2/1-49

AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA - Q9T0Q9_BPFD/1-49

FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA - COATB_BPIF1/22-73

你还可以使用比对信息对象的format 方法来以特定文件格式显示它–参见6.2.2 小节以了解细节。

你有没有注意到？在上面的原始文件中 ，很多序列 都包含咯数据库交叉引用信息，指向PDB 和关联的已知二级 （secondary）结构。试试看 ：

>>> for record in alignment:

...     if record.dbxrefs:

...         print record.id, record.dbxrefs

COATB_BPIKE/30-81 ['PDB; 1ifl ; 1-52;']

COATB_BPM13/24-72 ['PDB; 2cpb ; 1-49;', 'PDB; 2cps ; 1-49;']

Q9T0Q9_BPFD/1-49 ['PDB; 1nh4 A; 1-49;']

COATB_BPIF1/22-73 ['PDB; 1ifk ; 1-50;']

要看看所有序列的注释的话，试试这样：

>>> for record in alignment:

...     print record

Sanger在 http://pfam.sanger.ac.uk/family?acc=PF05371 提供咯一个漂亮的网页界面，它允许你以其它格式下载这个文件 。这是FASTA 格式的样子 ：

>COATB_BPIKE/30-81

AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIRLFKKFSSKA

>Q9T0Q8_BPIKE/1-52

AEPNAATNYATEAMDSLKTQAIDLISQTWPVVTTVVVAGLVIKLFKKFVSRA

>COATB_BPI22/32-83

DGTSTATSYATEAMNSLKTQATDLIDQTWPVVTSVAVAGLAIRLFKKFSSKA

>COATB_BPM13/24-72

AEGDDP---AKAAFNSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA

>COATB_BPZJ2/1-49

AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFASKA

>Q9T0Q9_BPFD/1-49

AEGDDP---AKAAFDSLQASATEYIGYAWAMVVVIVGATIGIKLFKKFTSKA

>COATB_BPIF1/22-73

FAADDATSQAKAAFDSLTAQATEMSGYAWALVVLVVGATVGIKLFKKFVSRA

注意，那个网页上应当有一个选项，控制的是把缺口 （gaps）显示成点还是减号，我们在上面是显示成减号的 。假设你下载咯它 ，并且保存成“PF05371_seed.faa”，那么你就可以使用几乎相同的代码来载入它咯 ：

from Bio import AlignIO

alignment = AlignIO.read("PF05371_seed.faa", "fasta")

print alignment

这个代码里面的变化仅仅是文件名和格式字符串。你会得到相同的输出 ，序列和记录标识符都是相同的。然而 ，你可能也会 猜到，如果你查看每个 SeqRecord 的话 ，将会发现没有注释也没有数据库交叉引用 ，因为FASTA 文件格式中没有这些东西。

注意，除咯使用Sanger 网页之外 ，你还可以使用 Bio.AlignIO 来亲自将原来的Stockholm 格式的文件转换成 FASTA 文件 （看下面的内容）。

对于任意受支持的文件格式，你都可以使用完全相同的方式来载入一条比对信息，仅仅需要改变格式字符串 。比如 ，使用“phylip”来读取PHYLIP 文件 、使用“nexus”来读取NEXUS 文件、使用“emboss”来读取EMBOSS 工具输出的比对文件。在维基页面 （ http://biopython.org/wiki/AlignIO ） 和内置文档 （也可 在线看 ）中有一个完整的列表：

>>> from Bio import AlignIO

>>> help(AlignIO)

...

6.1.2  多条比对信息

前一小节关注的是读取只包含单条比对信息的文件。然而 ，通常，文件中可以包含多条比对信息，要读取这些文件的话 ，我们必须使用 Bio.AlignIO.parse() 函数。

假设你有一条很短的比对信息，是PHYLIP 格式的：

    5    6

Alpha     AACAAC

Beta      AACCCC

Gamma     ACCAAC

Delta     CCACCA

Epsilon   CCAAAC

如果你想要使用PHYLIP 工具来引导（bootstrap）一个系统发生史 （phylogenetic） 树的话，那么 ， 其中的一个步骤就是 ，使用 bootseq 工具创建一组重采样 （resampled）过的比对信息。输出的内容是这样的 ，其中省略咯一些东西：

    5     6

Alpha     AAACCA

Beta      AAACCC

Gamma     ACCCCA

Delta     CCCAAC

Epsilon   CCCAAA

    5     6

Alpha     AAACAA

Beta      AAACCC

Gamma     ACCCAA

Delta     CCCACC

Epsilon   CCCAAA

    5     6

Alpha     AAAAAC

Beta      AAACCC

Gamma     AACAAC

Delta     CCCCCA

Epsilon   CCCAAC

...

    5     6

Alpha     AAAACC

Beta      ACCCCC

Gamma     AAAACC

Delta     CCCCAA

Epsilon   CAAACC

如果你想使用Bio.AlignIO 来读入这个的话，你可以这样：

from Bio import AlignIO

alignments = AlignIO.parse("resampled.phy", "phylip")

for alignment in alignments:

    print alignment

    print

将会产生以下输出，也是省略咯一些东西的：

SingleLetterAlphabet() alignment with 5 rows and 6 columns

AAACCA Alpha

AAACCC Beta

ACCCCA Gamma

CCCAAC Delta

CCCAAA Epsilon

SingleLetterAlphabet() alignment with 5 rows and 6 columns

AAACAA Alpha

AAACCC Beta

ACCCAA Gamma

CCCACC Delta

CCCAAA Epsilon

SingleLetterAlphabet() alignment with 5 rows and 6 columns

AAAAAC Alpha

AAACCC Beta

AACAAC Gamma

CCCCCA Delta

CCCAAC Epsilon

...

SingleLetterAlphabet() alignment with 5 rows and 6 columns

AAAACC Alpha

ACCCCC Beta

AAAACC Gamma

CCCCAA Delta

CAAACC Epsilon

就像使用Bio.SeqIO.parse()函数一样，使用 Bio.AlignIO.parse() 的时候也会返回一个迭代器。如果你想同时在内存中保存全部的比对信息，以 便以任意顺序访问它们的话 ，就把那个迭代器转换成一个列表：

from Bio import AlignIO

alignments = list(AlignIO.parse("resampled.phy", "phylip"))

last_align = alignments[-1]

first_align = alignments[0]

6.1.3  歧义性比对信息

很多比对文件格式都能显式储存多条比对信息，并且明显地标出各条比对信息之间的间隔。然而 ，如果使用的是一个通用的序列文件格式的话 ，就没有这种块状的结构 。最常见的情形就是使用FASTA 文件来存储比对信息的情况 。比如 ，看看下面的内容：

>Alpha

ACTACGACTAGCTCAG--G

>Beta

ACTACCGCTAGCTCAGAAG

>Gamma

ACTACGGCTAGCACAGAAG

>Alpha

ACTACGACTAGCTCAGG--

>Beta

ACTACCGCTAGCTCAGAAG

>Gamma

ACTACGGCTAGCACAGAAG

这个东西有可能是一个包含6个序列的单条比对信息（其中有重复的标识符）。或者 ，根据标识符判断，这可能是 2个不同的比对信息，每个都有3个序列，并且它们的长度碰巧相同。

下一个例子勒？

>Alpha

ACTACGACTAGCTCAG--G

>Beta

ACTACCGCTAGCTCAGAAG

>Alpha

ACTACGACTAGCTCAGG--

>Gamma

ACTACGGCTAGCACAGAAG

>Alpha

ACTACGACTAGCTCAGG--

>Delta

ACTACGGCTAGCACAGAAG

同样地，这个东西有可能是一个包含 6个序列的单条比对信息 。然而 ，根据标识符来看，我们可以猜测 ，这是 3个成对的比对信息，碰巧长度都一样。

最后勒个例子也类似：

>Alpha

ACTACGACTAGCTCAG--G

>XXX

ACTACCGCTAGCTCAGAAG

>Alpha

ACTACGACTAGCTCAGG

>YYY

ACTACGGCAAGCACAGG

>Alpha

--ACTACGAC--TAGCTCAGG

>ZZZ

GGACTACGACAATAGCTCAGG

在第三个例子中，由于长度不同，所以不能当成包含 6个序列的单条比对信息。然而 ，它可能是3条成对的比对信息。

显然，在一个FASTA 文件里保存多条比对信息是不合适的 。然而 ，如果你不得不使用这种输入文件 的话， Bio.AlignIO可能对付最普通的情况 ，那种情况下所有的比对信息 都有相同的记录个数。一个例子就是 一系列的成对比对信息 ，它们可能是由EMBOSS 工具needle 和water 生成的 –当然，在这种情况下， Bio.AlignIO应当能够理解它们的原配 （native）输出格式，使用“emboss”作为格式字符串。

要将这些FASTA 文件例子解释成多个单独的比对信息的话，我们可以使用 Bio.AlignIO.parse() ，并且使用可选的 seq_count 参数，它指定的是预期在每个比对信息中出现多少个序列 （在这些例子中 ，分别是 3 、 2和2）。比如 ，用第三个例子做输入数据 ：

for alignment in AlignIO.parse(handle, "fasta", seq_count=2):

    print "Alignment length %i" % alignment.get_alignment_length()

    for record in alignment:

        print "%s - %s" % (record.seq, record.id)

    print

输出：

Alignment length 19

ACTACGACTAGCTCAG--G - Alpha

ACTACCGCTAGCTCAGAAG - XXX

Alignment length 17

ACTACGACTAGCTCAGG - Alpha

ACTACGGCAAGCACAGG - YYY

Alignment length 21

--ACTACGAC--TAGCTCAGG - Alpha

GGACTACGACAATAGCTCAGG - ZZZ

对于前2个例子，使用 Bio.AlignIO.read() 或者不带 seq_count 参数的 Bio.AlignIO.parse() 的话 ，会输出单条比对信息 ，其中包含6个序列。对于第三个例子，会抛出一个异常，因为它们的长度不同 ，无法组成一个单条比对信息。

如果那个文件格式自身有一个块结构，使得 Bio.AlignIO 能够直接确定每个比对信息里的序列个数的话，那么就不需要 seq_count 参数。如果提供咯这个参数 ，却又与文件内容中的不一致的话，就会产生一个错误。

注意，这个可选的 seq_count 参数假设文件中的每个比对信息中有相同个数的序列 。假设 你可能遇到奇怪的情况，比如说一个包含 拥有不同序列个数的多条比对信息的FASTA 文件–当然，我很想能听说在真实世界中有这样的例子。假设你无法弄到一个有良好的文件格式的数据的话 ，那就没有简单的使用 Bio.AlignIO 来处理的方法 。在这种情况下 ，你可以考虑使用 Bio.SeqIO 将那些序列读入并且统一处理一下 ，再创建适当的比对文件 。