Biopython学习笔记

发表于 | 分类于

1. 序列读入与输出

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

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

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

  • Seq键：保存序列组成
  • ID键：保存序列ID

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

    ……

1.1. 解析/读取序列

  • （1）基本用法

    1
    SeqIO.parse( [FILE|HANDLE] , filetype)
    • 第一个参数是一个文件名或者一个句柄（ handle ，推荐）

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

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

      1
      2
      3
      4
      5
      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() 函数，如下：

    1
    2
    3
    4
    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 对象保持则一个列表，则可以这样：

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

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

  • （4）读取压缩文件里的序列

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

    1
    2
    3
    4
    5
    import 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
    2
    import 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()

    1
    orchid_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
    12
    # 定义key_function
    def get_accession(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]

    # 将SeqRecord对象转换为字典，并用前面定义的key_function定义字典的键
    orchid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.fasta", "fasta"), key_function=get_accession)

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

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

  • Bio.SeqIO.index()

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

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

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

1.3. 写入序列文件

Bio.SeqIO.write() 用于输出序列（写入文件）

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

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

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

1
2
3
4
from Bio import SeqIO
records = SeqIO.parse("ls_orchid.gbk", "genbank")
count = SeqIO.write(records, "my_example.fasta", "fasta")
print "Converted %i records" % count

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

1
count = 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
25
    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

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

然而，例如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
5
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() 写出序列比对文件

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

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

1
count = 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
10
from Bio import AlignIO
alignment = AlignIO.read("PF05371_seed.sth", "stockholm")
name_mapping = {}
for i, record in enumerate(alignment):
	# 建立一个字典对象实现自定义的ID和原始ID的映射
	name_mapping[i] = record.id
	record.id = "seq%i" % i
print name_mapping

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

2.3. 序列比对的操纵

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

1
alignment[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来执行序列比对任务的实现思路：

创造一个命令行对象来指定各种参数（例如：输入文件名，输出文件名等），然后通过Python的系统命令模块来运行这一程序

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

1
2
3
4
5
>>> import Bio.Align.Applications
>>> dir(Bio.Align.Applications)
...
['ClustalwCommandline', 'DialignCommandline', 'MafftCommandline', 'MuscleCommandline',
'PrankCommandline', 'ProbconsCommandline', 'TCoffeeCommandline' ...]

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

1
2
3
4
5
6
7
8
import os
# Python中 '\n' 和 '\t' 会被解析为一个新行和制表空白（tab），
# 如果你将一个小写的“r”放在字符串的前面，这一字符串就将保留原始状态，而不被解析
# 这种方式对于指定Windows风格的文件名来说是一种良好的习惯
clustalw_exe = r"C:\Program Files\new clustal\clustalw2.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

参考资料：

(1) BioPython中文官方文档