序列和序列对象

Seq 类

Seq类是Biopython最基础的一类, 储存序列信息. from Bio.Seq import Seq. 该类基本格式是Seq(self, data, alphabet=Alphabet()). 类似于字符串, 能够储存蛋白, DNA, RNA序列信息. 该类是不可变的. 该类和str类似, 支持count, find, split, strip.

相比str, Seq多了一个重要属性alphabet, 用于储存该序列的类型, 来源于Bio.Alphabet类似. 除此, 对于DNA序列还有如complement,

reverse_complement, transcribe, back_transcribe和translate 等特殊方法.

其基本使用例如:

>>> from Bio.Seq import Seq

>>> from Bio.Alphabet import IUPAC

>>> my_seq = Seq("MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF",

... IUPAC.protein)

>>> my_seq

Seq('MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF', IUPACProtein())

>>> print(my_seq)

MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVF

>>> my_seq.alphabet

IUPACProtein()

Seq类的特点如下:

具有alphabet属性, 储存了序列类型. 在切取序列后得到的新序列会保持该属性.

不同alphabet类型是不能互相+连接的. 如果要进行类似操作, 需要将alphabet更改为相同类型, 或者变为通用字母表.

具有字符串特性, 例如不可变, 能打印, 具有strip, split, find, startswith, endswith, lstrip等方法. 支持索引和切片操作.

count(substr)方法可以统计序列里字符或者片段的数量, 比较重要的方法.

lower(), upper() 方法可以变更大小写. IUPAC字母表均是大写.

str(Seq对象)或tostring()方法 可以转为str类型.

在进行in 判断, find方法等时, 大小写是敏感的.

alphabet 属性 (字母表)

字母表最重要的作用是标明序列里各个字母的含义, 标明序列的类型. 例如ATCG是DNA?RNA?还是蛋白呢?? 一般最好能够更为明确alphabet属性.

IUPAC字母表

IUPAC字母表是比较常用的字母表, 全部均是大写型. 如果对IUPAC字母表的序列进行lower, 类型会变成通用型.

from Bio.Alphabet import IUPAC，含有IUPAC的各种蛋白DNA和RNA的基本定义，对于识别序列内信息很重要。

Alphabet 类, 基本类, 但是字母类型不太明确.

IUPACProtein类 ，含有20基本氨基酸 (IUPAC.protein)

ExtendedIUPACProtein类，包含20种氨基酸外其他氨基酸元素(IUPAC.extended_protein)

IUPACUnambiguousDNA类，提供基本字母(IUPAC.unambiguous_dna)

IUPACAmbiguousDNA，提供每种可能下的歧义字母(IUPAC.ambiguous_dna)

ExtendedIUPACDNA类，提供修饰后的碱基(IUPAC.extended_dna)

IUPACUnambiguousRNA、IUPACAmbiguousRNA类。针对RNA。

>>> from Bio.Seq import Seq

>>> from Bio.Alphabet import IUPAC

>>> dna_seq = Seq("ACGT", IUPAC.unambiguous_dna)

>>> dna_seq

Seq('ACGT', IUPACUnambiguousDNA())

>>> dna_seq.lower()

Seq('acgt', DNAAlphabet())

通用字母表

Bio.Alphabet.generic_nucleotide, generic_dna, generic_rna, generic_protein, generic_alphabet, single_letter_alphabet.

这些是一些通用字母表. 通用字母表可以对下兼容, 例如generic_nucleotide类型序列可以和IUPAC.unambiguous_dna类型序列进行相加, 获得的是generic_nucleotide类型序列.

Seq序列的特殊方法

这些特殊方法可以在Bio.Seq内作为函数进行调用, 可以对字符串都进行相应处理. 例如:

>>> from Bio.Seq import reverse_complement, transcribe, back_transcribe, translate

>>> my_string = "GCTGTTATGGGTCGTTGGAAGGGTGGTCGTGCTGCTGGTTAG"

>>> reverse_complement(my_string)

'CTAACCAGCAGCACGACCACCCTTCCAACGACCCATAACAGC'

>>> transcribe(my_string)

'GCUGUUAUGGGUCGUUGGAAGGGUGGUCGUGCUGCUGGUUAG'

>>> back_transcribe(my_string)

'GCTGTTATGGGTCGTTGGAAGGGTGGTCGTGCTGCTGGTTAG'

>>> translate(my_string)

'AVMGRWKGGRAAG*'

核酸的complement(), 可以得到互补链.

核酸的reverse_complement() 可以得到反向的互补链.

DNA的transcribe()可以得到对应的mRNA的序列(其实就是T->U)

RNA的back_transcribe() 可以得到对应的DNA的序列(其实就是U->T)

核酸的translate() 可以将核酸翻译成蛋白序列, 包括DNA和RNA. 支持指定翻译表.

翻译结果里的*代表终止符.

coding_dna.translate(table="Vertebrate Mitochondrial") 可以指定使用线粒体的遗传密码表. 这个表也可以用table=2来指定. (参考NCBI基因编码)

to_stop=True形参可以使得遇到第一个终止密码子就停止, 并且不打印出来.

stop_symbol="@" 可以自己指定终止符.

cds=True可以说明该序列是完整CDS (首个会作为起始密码子, 可能翻译会不同)

>>> from Bio.Seq import Seq

>>> from Bio.Alphabet import IUPAC

### 互补链示例

>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)

>>> my_seq

Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())

>>> my_seq.complement()

Seq('CTAGCTACCCGGATATATCCTAGCTTTTAGCG', IUPACUnambiguousDNA())

>>> my_seq.reverse_complement()

Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC', IUPACUnambiguousDNA())

### 转录示例

>>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG", IUPAC.unambiguous_dna)

messenger_rna = coding_dna.transcribe()

>>> messenger_rna

Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG', IUPACUnambiguousRNA())

>>> messenger_rna.back_transcribe()

Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG', IUPACUnambiguousDNA())

### 翻译

>>> messenger_rna.translate()

Seq('MAIVMGR*KGAR*', HasStopCodon(IUPACProtein(), '*'))

>>> coding_dna.translate()

Seq('MAIVMGR*KGAR*', HasStopCodon(IUPACProtein(), '*'))

>>> coding_dna.translate(table="Vertebrate Mitochondrial")

Seq('MAIVMGRWKGAR*', HasStopCodon(IUPACProtein(), '*'))

翻译表

关于翻译表: 翻译表在CodonTable内, 其使用可以参考下面的例子:

>>> from Bio.Data import CodonTable

>>> standard_table = CodonTable.unambiguous_dna_by_name["Standard"]

>>> mito_table = CodonTable.unambiguous_dna_by_name["Vertebrate Mitochondrial"]

>>> print standard_table

........

>>> mito_table.stop_codons

['TAA', 'TAG', 'AGA', 'AGG']

>>> mito_table.start_codons

['ATT', 'ATC', 'ATA', 'ATG', 'GTG']

>>> mito_table.forward_table["ACG"]

'T'

序列的等同性

序列等同性比较奇怪, 教程里面认为不同的序列是不同的. 但实际测试, 即使是DNA和蛋白对比, 都会返回True (但会有Warning).

>>> from Bio.Seq import Seq

>>> from Bio.Alphabet import IUPAC

>>> seq1 = Seq("ACGT", IUPAC.ambiguous_dna)

>>> seq2 = Seq("ACGT", IUPAC.unambiguous_dna)

>>> seq3 = Seq("ACGT", IUPAC.protein)

>>> seq1 == seq2

True

>>> seq1 == seq3

# BiopythonWarning: Incompatible alphabets IUPACAmbiguousDNA() and IUPACProtein()

True

>>> str(seq1) == str(seq2)

True

MutableSeq 类型

Seq类型是不可变得, MutableSeq对象则是可变的. 其用法和Seq类似, 但可以改变内容. 该类型也可以由Seq类型的tomutable()方法进行转换.

MutableSeq 可变序列, 是不可以哈希化作为字典的键的, Seq则是可以.

该类型自然带一些新的方法:

append('A') : 可以在序列最后添加一个字符, 会改变原字符.

extend('TTTT') : 可以在序列最后添加新的内容. 可以延长一段序列.

insert(5, 'C') : 可以在相应索引位置插入内容.只能插入一个字符.

remove('A') : 可以在序列中移除首个遇到的指定字符.

pop() : 弹出序列最后一个字符(会改变原序列), 返回该字符.

reverse() : 反向链, 会修改原序列

reverse_complement() : 反向互补链, 会修改原序列.

>>> from Bio.Seq import Seq

>>> from Bio.Seq import MutableSeq

>>> from Bio.Alphabet import IUPAC

>>> my_seq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA", IUPAC.unambiguous_dna)

>>> mutable_seq = MutableSeq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA", IUPAC.unambiguous_dna)

>>> mutable_seq2 = my_seq.tomutable()

>>> mutable_seq2 MutableSeq('GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA', IUPACUnambiguousDNA())

>>> mutable_seq[5] = "C"

>>> mutable_seq.remove("T")

>>> mutable_seq.reverse()

UnknownSeq对象

UnknownSeq是Seq对象的一个子类, 目的是一个已知长度的序列, 但序列并不是由实际字母组成. 例如创造一个长度100万的"N"字母的序列, 但十分省内存.

>>> from Bio.Seq import UnknownSeq

>>> unk = UnknownSeq(20)

>>> unk

UnknownSeq(20, alphabet = Alphabet(), character = '?')

>>> print unk

????????????????????

>>> len(unk)

20

>>> unk_dna = UnknownSeq(20, alphabet=IUPAC.ambiguous_dna)

>>> unk_protein = unk_dna.translate()

>>> unk_protein

UnknownSeq(6, alphabet = ProteinAlphabet(), character = 'X')

SeqRecord 对象

from Bio.SeqRecord import SeqRecord. 一般使用SeqIO来读取创建, 也可以自行新建, SeqRecord 最少只需包含 Seq 对象.

.seq 序列自身(即 Seq 对象)。

.id 序列主ID(-字符串类型)。通常类同于accession number。

.name 序列名/id (-字符串类型)。 可以是accession number, 也可是clone名(类似GenBank record中的LOCUS id)。

.description : 序列描述(-字符串类型)。

.letter_annotations : 对照序列的每个字母逐字注释(per-letter-annotations)，以信息名为键(keys)，信息内容为值(value)所构成的字典。值与序列等长，用Python列表、元组或字符串表示。.letter_annotations可用于质量分数(如第 18.1.6 节) 或二级结构信息 (如 Stockholm/PFAM 比对文件)等数据的存储。

.annotations : 用于储存附加信息的字典。信息名为键(keys)，信息内容为值(value)。用于保存序列的零散信息(如unstructured information)。

.features : SeqFeature 对象列表，储存序列的结构化信息(structured information)，如：基因位置, 蛋白结构域。features 详见本章第三节( 第 4.3 节)。

.dbxrefs : 储存数据库交叉引用信息(cross-references)的字符串列表。

>>> from Bio import SeqIO

>>> record = SeqIO.read("NC_005816.fna", "fasta")

格式化方法

可以将记录转为SeqIO支持的格式, 如FASTA等. 例如record.format("fasta")

SeqFeature类对象

from Bio.SeqFeature import SeqFeature, FeatureLocation

SeqIO 类

Bio.SeqIO.parse()，它用于读取序列文件生成 SeqRecord对象，包含两个参数：

文件名或句柄.

小写字符串, 用于指定序列格式. 支持的格式参考Introduction to SeqIO. 读取的常见类型: "fasta", "genbank".

另外还有一个指定字母表的alphabet参数.

实际上, SeqIO.parse()函数返回时是SeqRecord对象迭代器, 可用于循环当中. 可用.next()方法遍历序列的条目 (最后没有条目时, 抛出StopIteration异常).

如果文件只有一个序列条目时, 可以使用SeqIO.read() , 具有类似参数, 返回一个SeqRecord对象而非迭代器.

>>> from Bio import SeqIO

>>> identifiers = [seq_record.id for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")]

>>> identifiers

['Z78533.1', 'Z78532.1', 'Z78531.1', 'Z78530.1', 'Z78529.1', 'Z78527.1', ..., 'Z78439.1']

压缩文件处理

如果下载的是压缩文件, 可以使用句柄的方法来读取文件:

>>> from Bio import SeqIO

## 对于gz压缩

>>> import gzip

>>> handle1 = gzip.open("ls_orchid.gbk.gz", "r")

>>> print sum(len(r) for r in SeqIO.parse(handle1, "gb"))

67518

>>> handle1.close()

## 对于bz2压缩

>>> import bz2

>>> handle2 = bz2.BZ2File("ls_orchid.gbk.bz2", "r")

>>> print sum(len(r) for r in SeqIO.parse(handle2, "gb"))

67518

>>> handle2.close()

网络序列条目

from Bio import Entrez

from Bio import SeqIO

Entrez.email = "A.N.Other@example.com"

handle = Entrez.efetch(db="nucleotide", rettype="fasta", retmode="text", id="6273291")

seq_record = SeqIO.read(handle, "fasta")

handle.close()

print "%s with %i features" % (seq_record.id, len(seq_record.features))

SwissPort序列条目

from Bio import ExPASy

from Bio import SeqIO

handle = ExPASy.get_sprot_raw("O23729")

seq_record = SeqIO.read(handle, "swiss")

handle.close()

print seq_record.id

print seq_record.name

print seq_record.description

print repr(seq_record.seq)

print "Length %i" % len(seq_record)

print seq_record.annotations["keywords"]

序列文件作为字典

SeqIO.to_dict(SeqRecord迭代器) : 获得一个字典, 默认键名是record的ID. 该方法适合较小的数据, 因为全部加载到内存! 如果想更改键名获取方法, 可以指定key_function来指定某个函数来获取键名.

SeqIO.index(文件名, 类型) : 类似于上一种,返回一个类似字典的对象, 但实际只记录每条序列条目在文件中的位置, 只在读取特定条目时才会进行解析. 输入是文件名, 类型, key_function等, 而不是SeqRecord迭代器. 返回的对象可以用get_raw()来提取原始数据.

SeqIO.index_db() : 将序列信息以文件方式存储到硬盘(SQLite3数据库), 因此适合处理超大文件. 该方法生成"索引文件idx", 返回的对象也支持get_raw()获取原始文件内容.

>>> from Bio import SeqIO

## todict方法

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

## index方法

>>> orchid_dict = SeqIO.index("ls_orchid.fasta", "fasta")

>>> len(orchid_dict)

94

>>> orchid_dict.keys()

['gi|2765596|emb|Z78471.1|PDZ78471', 'gi|2765646|emb|Z78521.1|CCZ78521', ...

..., 'gi|2765613|emb|Z78488.1|PTZ78488', 'gi|2765583|emb|Z78458.1|PHZ78458']

#### 使用index方法, 再用get_raw() 获取字符串格式的原始数据.

>>> from Bio import SeqIO

>>> uniprot = SeqIO.index("uniprot_sprot.dat", "swiss")

>>> handle = open("selected.dat", "w")

>>> for acc in ["P33487", "P19801", "P13689", "Q8JZQ5", "Q9TRC7"]:

... handle.write(uniprot.get_raw(acc))

>>> handle.close()

## index_db方法

>>> files = ["gbvrl%i.seq" % (i+1) for i in range(16)]

####### 保存到gbvrl.idx 索引文件(实际是SQLite3数据库).

>>> gb_vrl = SeqIO.index_db("gbvrl.idx", files, "genbank")

>>> print "%i sequences indexed" % len(gb_vrl)

958086 sequences indexed

>>> print gb_vrl["GQ333173.1"].description

HIV-1 isolate F12279A1 from Uganda gag protein (gag) gene, partial cds.

写入序列文件

使用Bio.SeqIO.write 来写出序列文件. 第一个参数是SeqRecord的列表或者迭代器(较新版本支持单个记录), 第二个参数是句柄或者输出文件名, 第三个参数是小写的输出格式.

使用Bio.SeqIO.convert()可以进行格式的转换, 而不用先打开再保存. 前两参数是输入文件名和格式, 后两参数是输出文件名和格式. 也可以用句柄. 注意, 由信息更少的格式转换为信息更多的格式时可能会缺失信息(例如FASTA->FASTAQ, 缺失质量值).

from Bio import SeqIO

## write的用法

my_records = [rec1, rec2, rec3]

SeqIO.write(my_records, "my_example.faa", "fasta")

## convert转换格式的用法

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

print "Converted %i records" % count

NCBI运行BLAST

使用Bio.Blast.NCBIWWW 模块的qblast()来调用在线版本的BLAST, 最后返回一个句柄. 这个函数含有3个必要的参数:

搜索的blast程序, 小写字符串.

指定搜索的数据库, 可参考 搜索数据库

查询的字符串, 可以是序列本身(fasta格式)或者类似GI的id.

支持的blast程序

程序

功能介绍

blastn

搜索核酸 : 在核酸库

blastp

搜索蛋白: 在蛋白库

blastx

搜索核酸: 在蛋白库, 在所有6种框架下动态翻译

tblastn

搜索蛋白: 在核酸库, 在所有6种框架下动态翻译

tblastx

搜索核酸: 在翻译的核酸库, 在所有6种框架下动态翻译

除此以外, 还支持其他一些参数:

format_type可以指定返回格式关键字, 如"HTML", "Text", "ASN.1", 或 "XML". 默认是XML.

expect, 即阈值 e-value.

## 通过GI来指定

>>> from Bio.Blast import NCBIWWW

>>> result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")

## 通过SeqRecord 来指定

>>> from Bio import SeqIO

>>> record = SeqIO.read("m_cold.fasta", format="fasta")

>>> result_handle = NCBIWWW.qblast("blastn", "nt", record.format("fasta"))

本地运行BLAST

使用Bio.Blast.Applications的NcbiblastxCommandline来构建命令行字符串并运行它(函数运行). 另外, 还有blastx替换为blastn, blastp, deltablast, psiblast,rpsblast, rpstblastn, tblastn,tblastx, blastformatter的相应程序.

>>> from Bio.Blast.Applications import NcbiblastxCommandline

>>> help(NcbiblastxCommandline)

...

>>> blastx_cline = NcbiblastxCommandline(query="opuntia.fasta", db="nr", evalue=0.001, outfmt=5, out="opuntia.xml")

>>> blastx_cline

NcbiblastxCommandline(cmd='blastx', out='opuntia.xml', outfmt=5, query='opuntia.fasta', db='nr', evalue=0.001)

>>> print blastx_cline

blastx -out opuntia.xml -outfmt 5 -query opuntia.fasta -db nr -evalue 0.001

>>> stdout, stderr = blastx_cline()

## 注意这里用了输出到XML, 因此stdout和stderr是空的.

解析BLAST输出

最好用也最推荐的输出格式是XML. 无论是网页运行, 本地运行, 还是Biopython运行, 都可以指定生成该格式, 最好使用句柄打开和处理他. 类似之前SeqIO或AlignIO, 有read和parse两个方法处理结果, 前者针对一个结果, 后者针对多个结果.

### 处理器核心 NCBIXML

>>> from Bio.Blast import NCBIXML

### 获取一份结果

>>> from Bio.Blast import NCBIWWW

>>> result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")

#>>> result_handle = open("my_blast.xml")

### 处理结果

>>> blast_record = NCBIXML.read(result_handle)

### 如果用了多条序列去blast, 有许多的搜索结果:

>>> blast_records = NCBIXML.parse(result_handle)

使用read或者parse解析后获得的对象是Bio.Blast.Record.Blast 类对象. 内有丰富的信息.

.alignments : 重要的aligment结果的list. 找到的结果都在这.

.descriptions : 描述

.database, expect : 搜索的数据库, e值,

.query, .query_id, .query_length 搜索的序列信息, ID, 长度.

重点解析的内容是alignments的结果. alignment含有的属性有:

.hsps : HSP类对象的列表list. 结果主要信息.

.accession : Accession 号, gi|1219041180|ref|XM_021875076.1|后面的XM_021875076.

.length : 该序列的长度.

.title : 该序列的头信息.

.hit_def : 该title最后的叙述信息.

.hit_id : 该title的前面ID信息.

对于 HSP类对象( high-scoring pair(高分片段), 含有以下属性:

.score, bits, .expect : 实际的blast分值, 分值相应的匹配数, e值.

.query : 查询的序列的匹配部分, query_start/end, 匹配部分起始号.

.match : 匹配的情况, |的是匹配得上的部分

.sbjct : 匹配的序列的匹配部分, sbjct_start/end, 匹配部分起始号.

.gaps : gaps总数 (针对query)

.identities, .positives : 匹配部分相同/正分的数量.

>>> E_VALUE_THRESH = 0.04

>>> for alignment in blast_record.alignments:

... for hsp in alignment.hsps:

... if hsp.expect < E_VALUE_THRESH:

... print '****Alignment****'

... print 'sequence:', alignment.title

... print 'length:', alignment.length

... print 'e value:', hsp.expect

... print hsp.query[0:75] + '...'

... print hsp.match[0:75] + '...'

... print hsp.sbjct[0:75] + '...'

Bio.SeqUtils.GC(序列) 可以返回序列里GC的百分比。

Reference