怎么用python实现序列比对_第6章 多序列比对

6.1.1 单一的序列比对¶

例如，请见以下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

//

这是PFAM数据库中Phage_Coat_Gp8的种子排列(PF05371)。该排列下载于一个已经过期的PFAM数据库版本( http://pfam.sanger.ac.uk/ )，但这并不影响我们的例子(假设你已经将以上内容下载到一个名为’‘PF05371_seed.sth’‘的文件中，并在Python的当前工作目录下)：

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

你也可以使用上面alignment对象的 format 方法来以指定的格式显示它。具体信息可以参见 6.2.2 。

你是否已经注意到以上原始数据文件中包含有蛋白数据库(PDB)交叉引用以及相关二级结构的信息？你可以尝试以下代码：

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

>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

注意Sanger网站有一个选项可以将序列比对中的间隔(gap)用小圆点或者是小横线表示。在以上例子中，序列间隔由小横线表示。假设你已经下载该文件，并保存为 “PF05371_seed.faa”。你可以使用以下代码来读入该序列比对。

from Bio import AlignIO

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

print alignment

你可能已经发现，以上代码中唯一的变化只是指定格式的参数。所返回的alignment对象将会包含同样的序列和序列名字。但是仔细的读者会发现，每一个alignment的SeqRecord中并不包含数据的引用注释。这是因为FASTA格式本身并没有包含这一类信息。

此外，除了使用Sanger网站，你也可以利用 Bio.AlignIO 来将原始的Stockholm格式转换成FASTA文件格式(见下文)。

对于任何一种Biopython支持的格式，你都可以用同样的方式读取它(通过指定文件的格式)。例如，你可以使用“phylip”来表示PHYLIP格式文件，用”nexus”来指定NEXUS格式文件或者用“emboss”来指定EMBOSS工具箱的输出文件。读者可以在以下链接中找到所有支持的格式( http://biopython.org/wiki/AlignIO )，或者内置的帮助中(以及在线文档 online )：

>>>from Bio import AlignIO

>>>help(AlignIO)

...