biopython中的代码 biopython教程与手册

转载

小咪咪 2024-01-24 15:35:20

文章标签 biopython中的代码 biopython中文指南句柄迭代器自定义 文章分类 Python 后端开发

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 )，迭代器通常用在循环中

from 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() 函数，如下：

from 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 模块打开压缩文档以读取数据

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

mport 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 ，它允许你根据你的序列条目特点，自定义字典键，例如：

2def (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 迭代器 - 这允许我们通过结合使用这两个函数实现文档转换。

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格式的一个多序列比对文档：

# 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格式：

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一类的普通序列文档格式并没有很直接的分隔符来分开多个序列比对：

>Alpha
ACTACGACTAGCTCAG--G
>Beta
ACTACCGCTAGCTCAGAAG
>Gamma
ACTACGGCTAGCACAGAAG
>Alpha
ACTACGACTAGCTCAGG--
>Beta
ACTACCGCTAGCTCAGAAG
>Gamma
ACTACGGCTAGCACAGAAG

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

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

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

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

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的映射)：

from 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条序列

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

>>> print alignment[2,6]
T
# 其实是以下操作的简化
>>> print alignment[2].seq[6]
T
# 获取整列
>>> print alignment[:,6]
TTT---T

2.4. 构建序列比对的工具

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

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

import Bio.Align.Applications
dir(Bio.Align.Applications)
...
['ClustalwCommandline', 'DialignCommandline', 'MafftCommandline', 'MuscleCommandline',
'PrankCommandline', 'ProbconsCommandline', 'TCoffeeCommandline' ...]

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

import 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

>>> from Bio.Align.Applications import ClustalwCommandline
>>> cline = ClustalwCommandline("clustalw2", infile="opuntia.fasta")
>>> print cline
clustalw2 -infile=opuntia.fasta

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。