biopython中文指南_Biopython学习笔记

1. 序列读入与输出

序列的读入与输出通过Bio.SeqIO模块实现，它旨在提供一个简单的接口，实现对各种不同格式序列文档进行统一的处理

主要是基于对SeqRecord对象的操作，该对象包含一个 Seq 对象)和注释信息(如序列ID和描述信息)

SeqRecord对象本质上就是一个字典：

Seq键：保存序列组成

ID键：保存序列ID

description键：保存序列的描述信息

……

1.1. 解析/读取序列(1)基本用法

1SeqIO.parse( [FILE|HANDLE] , filetype)

第一个参数是一个文档名或者一个句柄( handle ，推荐)

第二个参数是一个小写字母字符串，用于指定序列格式(我们并不推测文档格式！)，支持的文档格式

Bio.SeqIO.parse() 函数返回一个 SeqRecord 对象迭代器( iterator )，迭代器通常用在循环中

1

2

3

4

5from Bio import SeqIO

for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):

print seq_record.id

print seq_record.seq

print len(seq_record)

(2)使用with和句柄的方法操作

由于我们建议养成一个良好的编程习惯——在文档句柄用完后要及时把它关闭，所以推荐使用with和句柄的方法操作Bio.SeqIO.parse() 函数，如下：

1

2

3

4from Bio import SeqIO

with open("ls_orchid.gbk") as handle:

for r in SeqIO.parse(handle, "gb"):

...

有时你需要处理只包含一个序列条目的文档，此时请使用函数 Bio.SeqIO.read() 。它使用与函数 Bio.SeqIO.parse() 相同的参数，因为Bio.SeqIO.parse() 函数返回一个 SeqRecord 对象迭代器，所以也可以使用Bio.SeqIO.parse().next()

(3)保存成列表形式——方便后续操作

如果想要把整个文档中的所有序列的SeqRecord 对象保持则一个列表，则可以这样：

1records = [seq_record for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")]

将序列保持成一个列表，使我们能方便地从里面挑取里面的任意一条序列的信息

(4)读取压缩文档里的序列

可以使用Python的 gzip 模块打开压缩文档以读取数据

1

2

3

4

5import gzip

from Bio import SeqIO

with gzip.open("ls_orchid.gbk.gz", "r") as handle:

for record in SeqIO.parse(handle,'gb'):

...

相同地，如果我们有一个bzip2压缩文档，可以使用bz2模块

1

2import bz2

handle = bz2.BZ2File("ls_orchid.gbk.bz2", "r")

1.2. 转换成字典形式

之所以要将序列文档的解析结果转换为字典形式，是为了后续能随意地读取任意一条序列的信息

Bio.SeqIO 模块中3个相关函数，用于随机读取多序列文档

Bio.SeqIO.to_dict()：最灵活但内存占用最大，因为它将整个字典都载人到内存中；

Bio.SeqIO.index()：工作原理上与Bio.SeqIO.to_dict()略有不同——它是对这个文档进行索引。尽管仍然是返回一个类似于字典的对象，它并不将所有的信息存储在内存中。相反，它仅仅记录每条序列条目在文档中的位置——当你需要读取某条特定序列条目时，它才进行解析；

Bio.SeqIO.to_dict()

1orchid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.gbk", "genbank"))

使用 Bio.SeqIO.to_dict() 函数将明确检查重复键，如果发现任何重复键将引发异常并退出

默认会使用每条序列条目的ID(i.e. .id 属性)作为键

如果你需要别的作为键，如登录号(Accession Number)，可使用 SeqIO.to_dict() 的可选参数 key_function ，它允许你根据你的序列条目特点，自定义字典键，例如：

1

2

3

4

5

6

7

8

9

10

11

12def (record):

""""Given a SeqRecord, return the accession number as a string.

e.g. "gi|2765613|emb|Z78488.1|PTZ78488" -> "Z78488.1"

"""

parts = record.id.split("|")

assert len(parts) == 5 and parts[0] == "gi" and parts[2] == "emb"

return parts[3]

orchid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.fasta", "fasta"), key_function=get_accession)

也可以利用python的lambda表达式定义一个简易的一次性的函数来充当key_function：

1seguid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.gbk", "genbank"), lambda rec : seguid(rec.seq))

Bio.SeqIO.index()

1orchid_dict = SeqIO.index("ls_orchid.gbk", "genbank")

注意： Bio.SeqIO.index() 不接受句柄参数，仅仅接受文档名

它默认也会使用每条序列条目的ID(i.e. .id 属性)作为键，如果想要自定义，也需要定义一个key_function

1.3. 写入序列文档

Bio.SeqIO.write() 用于输出序列(写入文档)

1Bio.SeqIO.write(SeqRecord, [FILE | HANDLE], file-format)

可以在写出时实现文档格式的转换：

在之前的例子中我们使用 SeqRecord 对象列表作为 Bio.SeqIO.write() 函数的输入，但是它也接受如来自于 Bio.SeqIO.parse() 的 SeqRecord 迭代器 - 这允许我们通过结合使用这两个函数实现文档转换。

1

2

3

4from Bio import SeqIO

records = SeqIO.parse("ls_orchid.gbk", "genbank")

count = SeqIO.write(records, "my_example.fasta", "fasta")

print "Converted %i records" % count

这仍然有点复杂，有一个辅助函数可以替代上述代码：

1count = SeqIO.convert("ls_orchid.gbk", "genbank", "my_example.fasta", "fasta")

2. 序列比对

2.1. 读取多序列比对数据

在Biopython中，有两种方法读取多序列比对数据：

Bio.AlignIO.read() 只能读取一个多序列比对；

Bio.AlignIO.parse() 可以依次读取多个序列比对数据。

使用 Bio.AlignIO.parse() 将会返回一个 MultipleSeqAlignment 的 迭代器(iterator)

Stockholm格式的一个多序列比对文档：

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

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

//

phylip格式：

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

255 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

这些格式的比对文档都可以非常明确地储存多个序列比对

然而，例如FASTA一类的普通序列文档格式并没有很直接的分隔符来分开多个序列比对：

1

2

3

4

5

6

7

8

9

10

11

12>Alpha

ACTACGACTAGCTCAG--G

>Beta

ACTACCGCTAGCTCAGAAG

>Gamma

ACTACGGCTAGCACAGAAG

>Alpha

ACTACGACTAGCTCAGG--

>Beta

ACTACCGCTAGCTCAGAAG

>Gamma

ACTACGGCTAGCACAGAAG

以上FASTA格式文档可以认为是一个包含有6条序列的序列比对(有重复序列名)。或者从文档名来看，这很可能是两个序列比对，每一个包含有三个序列，只是这两个序列比对恰好具有相同的长度

很明显，将多个序列比对以FASTA格式储存并不方便

为了处理这样的FASTA格式的数据，我们可以指定 Bio.AlignIO.parse() 的第三个可选参数 seq_count ，这一参数将告诉Biopython你所期望的每个序列比对中序列的个数。例如：

1

2

3

4

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

2.2. 序列比对的写出

使用 Bio.AlignIO.write() 写出序列比对文档

1Bio.AlignIO.write(MultipleSeqAlignment, [FILE | HANDLE], file-format)

Bio.AlignIO 模块中的序列比对格式转换功能与 Bio.SeqIO 模块的格式转换是一样的

1count = AlignIO.convert("PF05371_seed.sth", "stockholm", "PF05371_seed.aln", "clustal")

当你写出的文档格式是PHYLIP时，需要注意：

PHYLIP格式它严格地要求每一条序列的ID是都为10个字符(ID中多出的字符将被截短)，这会带来序列ID的改变，甚至使得原先ID不同的两条序列，被截短后ID变成一样的了

为了解决这个问题，我们提供了另一种解决方案 —— 利用自定义的序列ID来代替原本的序列ID(建立一个字典对象实现自定义的ID和原始ID的映射)：

1

2

3

4

5

6

7

8

9

10from Bio import AlignIO

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

name_mapping = {}

for i, record in enumerate(alignment):

name_mapping[i] = record.id

record.id = "seq%i" % i

print name_mapping

AlignIO.write([alignment], "PF05371_seed.phy", "phylip")

2.3. 序列比对的操纵

使用切片操作来获得其中某些序列比对

1alignment[3:7] # 获取一个多序列比对的第4到第8条序列

获得特定的列，使用双切片：

1

2

3

4

5

6

7

8

9>>> print alignment[2,6]

T

# 其实是以下操作的简化

>>> print alignment[2].seq[6]

T

# 获取整列

>>> print alignment[:,6]

TTT---T

2.4. 构建序列比对的工具

使用Python来执行序列比对任务的实现思路：

大多数的打包进程都在 Bio.Align.Applications 中定义：

1

2

3

4

5import Bio.Align.Applications

dir(Bio.Align.Applications)

...

['ClustalwCommandline', 'DialignCommandline', 'MafftCommandline', 'MuscleCommandline',

'PrankCommandline', 'ProbconsCommandline', 'TCoffeeCommandline' ...]

注意：在python调用这些序列比对工具之前，需要保证这些工具已经安装好，且被添加到PATH下，或者也可以提供这些工具的绝对路径，下面以ClustalW为例进行说明：

1

2

3

4

5

6

7

8import os

# 如果你将一个小写的“r”放在字符串的前面，这一字符串就将保留原始状态，而不被解析

# 这种方式对于指定Windows风格的文档名来说是一种良好的习惯

clustalw_exe = r"C:Program Filesnew clustalclustalw2.exe"

# 检测指定的路径是否存在

assert os.path.isfile(clustalw_exe), "Clustal W executable missing"

clustalw_cline = ClustalwCommandline(clustalw_exe, infile="opuntia.fasta")

ClustalW

1

2

3

4>>> from Bio.Align.Applications import ClustalwCommandline

>>> cline = ClustalwCommandline("clustalw2", infile="opuntia.fasta")

>>> print cline

clustalw2 -infile=opuntia.fasta

参考资料：