biopython

转载
Part 2  Biopython的重头戏-生物学中序列的处理

Biopyhton的Seq和Python中标准字符串有两大重要的不同之处：首先，他们的处理方法不同。Seq适用于很多不同字符串的用的方法，如translate（）,但是又有所不同。而且Biopython中加入了不同字符处理中没有的方法，如reverse_complement（）; 第二，Seq模块中加入了重要的特性alphabet，这个对象可以解释序列代表的意思，即这个序列是一个DNA序列还是蛋白质序等。
 是alphabet对象使得Seq对象成为不同于普通字符串的重要标志。Alphabet在Bio.Alphabet中定义。我们用IUPAC alphabet来处理很对重要的对象：DNA，RNA， Proteins.
Bio.Alphabet.IUPAC提供了写了的基本信息。比如对于常见蛋白质序列用IUPACProtein来定义，而对于人体中不常见的氨基酸序列则用ExtendedIUPACProtein来定义；DNA序列用IUPACUnambigousDNA来定义，对于特殊的DNA序列（如甲基化的DNA序列）则用ExtendedIUPACDNA来定义等等。
现在我们可以来处理生物学序列了，但是你要先装上Biopython模块。
 
1.   基本生物学序列的表示与处理
 
a)    对于模糊不清的序列我们可以这样处理
>>> from Bio.Seq import Seq  #import Seq 模块来处理序列
>>> my_seq = Seq("AGTACACTGGT")  #定义自己的序列Seq()
>>> my_seq
Seq(’AGTACACTGGT’, Alphabet())  #这里不是像string那样单单显示‘’AGTA..‘’了
#可以看到，这里的seq是包含Alphabet特性的，和普通序列不同
>>> my_seq.alphabet Alphabet()
my_seq.alphabet Alphabet()
 
b)    对于要详细说明序列，就要用到alphabet中的详细特性了。定义DNA序列：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC     #前面用介绍
>>> my_seq = Seq("AGTACACTGGT", IUPAC.unambiguous_dna)
>>> my_seq
Seq(’AGTACACTGGT’, IUPACUnambiguousDNA())
>>> my_seq.alphabet    #相当于显示序列属性
IUPACUnambiguousDNA()  
 
c)    定义蛋白质序列
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_prot = Seq("AGTACACTGGT", IUPAC.protein)
>>> my_prot
Seq(’AGTACACTGGT’, IUPACProtein())
>>> my_prot.alphabet
IUPACProtein()
 

2.   像处理普通字符串那样处理生物学的序列
 
a)    适用于普通String的方法，如
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("GATCG", IUPAC.unambiguous_dna)
>>> for index, letter in enumerate(my_seq):  #索引，print等特性
...print index, letter
0 G
1 A
2 T
3 C
4 G
>>> print len(my_seq)  #拥有像不同字符串的处理方法，很酷吧！
5
>>> print my_seq[0] #first letter
G
>>> print my_seq[2] #third letter
T
>>> print my_seq[-1] #last letter
G
>>> "AAAA".count("AA")
2
>>> Seq("AAAA").count("AA")
2
 
b）对于特殊的生物学问题处理起来就很方便，如计算GC%：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq(’GATCGATGGGCCTATATAGGATCGAAAATCGC’, IUPAC.unambiguous_dna)
>>> len(my_seq)
32
>>> my_seq.count("G")
9
>>> 100 * float(my_seq.count("G") + my_seq.count("C")) / len(my_seq)
46.875
 
C) 当然你可以用上面的方法计算GC%，但是Biopython中还有更方便的模块Bio.SeqUtils计算GC%
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> from Bio.SeqUtils import GC
>>> my_seq = Seq(’GATCGATGGGCCTATATAGGATCGAAAATCGC’,IUPAC.unambiguous_dna)
>>> GC(my_seq)   #很酷
46.875
 

3.   序列的分片
 
从序列中截取想要的序列就要分片了，如
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC",IUPAC.unambiguous_dna)
>>> my_seq[4:12]
Seq(’GATGGGCC’, IUPACUnambiguousDNA())  #得到的仍然是具有属性的
>>> my_seq[0::3]
Seq(’GCTGTAGTAAG’, IUPACUnambiguousDNA())  #从0开始，中间间隔3个去，下面同理
>>> my_seq[1::3]
Seq(’AGGCATGCATC’, IUPACUnambiguousDNA())
>>> my_seq[2::3]
Seq(’TAGCTAAGAC’, IUPACUnambiguousDNA())
 
4.   将生物学序列转化成普通字符串
 
方法1
>>> str(my_seq)
’GATCGATGGGCCTATATAGGATCGAAAATCGC’
>>> print my_seq
GATCGATGGGCCTATATAGGATCGAAAATCGC
 
方法2：将FASTA格式中的序列提取出来：
>>> fasta_format_string = ">Name\n%s\n" % my_seq
>>> print fasta_format_string
>Name
GATCGATGGGCCTATATAGGATCGAAAATCGC
 方法3：
>>> my_seq.tostring()
’GATCGATGGGCCTATATAGGATCGAAAATCGC’
 
5.   序列的拼接：
 
>>> from Bio.Alphabet import IUPAC
>>> from Bio.Seq import Seq
>>> protein_seq = Seq("EVRNAK", IUPAC.protein)
>>> dna_seq = Seq("ACGT", IUPAC.unambiguous_dna)
>>> protein_seq + dna_seq
Traceback (most recent call last):
TypeError: Incompatible alphabets IUPACProtein() and IUPACUnambiguousDNA()
#序列拼接的简单方法就是+，但是不同性质的序列是不可以相加的，所以DNA序列和蛋白序列相加会报错！可以像下面一样统一特性。
 
>>> from Bio.Alphabet import generic_alphabet
>>> protein_seq.alphabet = generic_alphabet
>>> dna_seq.alphabet = generic_alphabet
>>> protein_seq + dna_seq
Seq(’EVRNAKACGT’, Alphabet())   #这里的Alphabet是特性变得模糊了
 
另外对于一般的核苷酸序列特性和IUPAC DNA sequence相加也会使得结果变得模糊（为一般核苷酸序列特性）！如：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import generic_nucleotide
>>> from Bio.Alphabet import IUPAC
>>> nuc_seq = Seq("GATCGATGC", generic_nucleotide)
>>> dna_seq = Seq("ACGT", IUPAC.unambiguous_dna)
>>> nuc_seq
Seq(’GATCGATGC’, NucleotideAlphabet())
>>> dna_seq
Seq(’ACGT’, IUPACUnambiguousDNA())
>>> nuc_seq + dna_seq
Seq(’GATCGATGCACGT’, NucleotideAlphabet())
 
6.   改变大小写等情况
 
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import generic_dna
>>> dna_seq = Seq("acgtACGT", generic_dna)
>>> dna_seq
Seq(’acgtACGT’, DNAAlphabet())
>>> dna_seq.upper()
Seq(’ACGTACGT’, DNAAlphabet())
>>> dna_seq.lower()
Seq(’acgtacgt’, DNAAlphabet())
>>> "GTAC" in dna_seq
False
>>> "GTAC" in dna_seq.upper()
True
 
注意IUPAC alphabets仅仅识别大写字母：
>>> 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())
 
7.   核苷酸序列的反向，互补等方法
 
对于核苷酸序列，可以直接得到其反向互补的序列的。
>>> 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())
 
另外最普通的反转是这样：
>>> my_seq[::-1]
Seq(’CGCTAAAAGCTAGGATATATCCGGGTAGCTAG’, IUPACUnambiguousDNA())
 
蛋白序列没有反转互补等特性：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> protein_seq = Seq("EVRNAK", IUPAC.protein)
>>> protein_seq.complement()
Traceback (most recent call last):
...
ValueError: Proteins do not have complements!
 
 8.   重要的生物学过程之转录（DNA-mRNA）
 
在讨论转录之前，先让我们熟悉一下转录的基本知识。首先要分清编码链和模板链。mRNA是与模板链互补配对的，mRNA的序列是和编码链相同的；DNA中的T在RNA中相应的为U，即DNA中A的互补配对碱基是RNA中的U；一般情况下默认是通过模板链3‘-5’，转录成mRNA是5‘-3’。
 
以上可以总结为下图
DNA coding strand (aka Crick strand, strand +1)
5’ ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG 3’
3’ TACCGGTAACATTACCCGGCGACTTTCCCACGGGCTATC 5’
DNA template strand (aka Watson strand, strand −1)
 
Transcription
 
5’ AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG 3’
Single stranded messenger RNA
 
Biopython对于处理转录有很大的优势，如：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG", IUPAC.unambiguous_dna)
>>> coding_dna
Seq(’ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG’, IUPACUnambiguousDNA())
>>> template_dna = coding_dna.reverse_complement()    #编码链和模板链是互补配对的
>>> template_dna
Seq(’CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT’, IUPACUnambiguousDNA())
#应该是默认模板链的书写方式是从3‘-5’的。
 
我们知道mRNA中没有T，通过模板链直接reverse.complement()是不能得到结果的，这是就要用到seq.transcribe()方法了。
>>> messenger_rna = coding_dna.transcribe()
>>> messenger_rna
Seq(’AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG’, IUPACUnambiguousRNA())
 
还可以直接用一句话搞定：
>>> template_dna.reverse_complement().transcribe()
Seq(’AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG’,
IUPACUnambiguousRNA())
 
对于反转录，也有相应的方法，反转录后得到的编码链的序列。
>>> messenger_rna.back_transcribe()
Seq(’ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG’, IUPACUnambiguousDNA())
 
9.   重要的生物学过程之二翻译
 
翻译是从mRNA到proteins的过程。真核起始密码子AUG，终止密码子UAG，UGA，UAA。原核和真核是不同的。
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> messenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG", IUPAC.unambiguous_rna)
>>> messenger_rna
Seq(’AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG’, IUPACUnambiguousRNA())
>>> messenger_rna.translate()
Seq(’MAIVMGR*KGAR*’, HasStopCodon(IUPACProtein(), ’*’)) #*表示终止密码子，上个序列中存在两个终止密码子，不编码氨基酸。
 
也可以直接从DNA翻译成为氨基酸序列
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG", IUPAC.unambiguous_dna)
>>> coding_dna
Seq(’ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG’, IUPACUnambiguousDNA())
>>> coding_dna.translate()
Seq(’MAIVMGR*KGAR*’, HasStopCodon(IUPACProtein(), ’*’))
 
以上的翻译是安装NCBI中的翻译表进行转换的，当然你可以指定用其他的译码方式。但是Biopython默认的是standard genetic code (NCBI table id 1).你可以自己指定
如当你要处理线粒体的DNA时候，要这样：
>>> coding_dna.translate(table="Vertebrate Mitochondrial")  #当是细菌的时候，应为Bacterial
Seq(’MAIVMGRWKGAR*’, HasStopCodon(IUPACProtein(), ’*’)) #和默认的不同吧
当你想指定用另一个表时，可以这样设置：
>>> coding_dna.translate(table=2)
Seq(’MAIVMGRWKGAR*’, HasStopCodon(IUPACProtein(), ’*’))
 
当你想使遇到第一个密码子就停止翻译时（自然就是这样），你可以这样设置：
>>> coding_dna.translate()
Seq(’MAIVMGR*KGAR*’, HasStopCodon(IUPACProtein(), ’*’))
>>> coding_dna.translate(to_stop=True)
Seq(’MAIVMGR’, IUPACProtein())  #和前者不同吧
>>> coding_dna.translate(table=2)
Seq(’MAIVMGRWKGAR*’, HasStopCodon(IUPACProtein(), ’*’))
>>> coding_dna.translate(table=2, to_stop=True) #以表2翻译，遇到第一个终止密码子即终止
Seq(’MAIVMGRWKGAR’, IUPACProtein())
当然你也可以指定终止符号：
>>> coding_dna.translate(table=2, stop_symbol="@")
Seq(’MAIVMGRWKGAR@’, HasStopCodon(IUPACProtein(), ’@’))
 
然而当你的转录起始序列不是AUG时，即其他生物的密码和真核可能不一样，这时候也要设置：
>>> gene.translate(table="Bacterial", to_stop=True)
 
10. 翻译表
 
这里我们紧紧关注两类翻译表：标准翻译表和脊椎动物线粒体转录表
>>> from Bio.Data import CodonTable
>>> standard_table = CodonTable.unambiguous_dna_by_name["Standard"]
>>> mito_table = CodonTable.unambiguous_dna_by_name["Vertebrate Mitochondrial"]
上面的命令代码和下面是等价的
>>> from Bio.Data import CodonTable
>>> standard_table = CodonTable.unambiguous_dna_by_id[1]
>>> mito_table = CodonTable.unambiguous_dna_by_id[2]
 
你可以将表显示出来：
>>> 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’
 

11.生物序列的比较：
 
这里的序列和普通的序列不同，不尽要考虑序列还要考虑Alphabet特性
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> seq1 = Seq("ACGT", IUPAC.unambiguous_dna)
>>> seq2 = Seq("ACGT", IUPAC.unambiguous_dna)
>>> seq1 == seq2  #这个为什么不对呢，我也不知道
False
>>> seq1 == seq1
True
>>> id(seq1) == id(seq2)
False
>>> id(seq1) == id(seq1)
True
>>> str(seq1) == str(seq2)
True
>>> str(seq1) == str(seq1)
True
 
12.可变的对象：
 
当你想改变一个序列中的序列时，怎么办呢
你可以先试试：
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA", IUPAC.unambiguous_dna)
>>> my_seq[5] = "G"
Traceback (most recent call last):
...
TypeError: ’Seq’ object does not support item assignment     
#报错了，这种直接的方法是不可行的，试试下面的方法：
>>> mutable_seq = my_seq.tomutable()   #使得序列可变的一种方法。
>>> mutable_seq
MutableSeq(’GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA’, IUPACUnambiguousDNA())
Alternatively, you
 
当然你也可用下边一种方法直接生成可变序列：
>>> from Bio.Seq import MutableSeq
>>> from Bio.Alphabet import IUPAC
>>> mutable_seq = MutableSeq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA", IUPAC.unambiguous_dna)
>>> mutable_seq
MutableSeq(’GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA’, IUPACUnambiguousDNA())
>>> mutable_seq[5] = "C"
>>> mutable_seq
MutableSeq(’GCCATCGTAATGGGCCGCTGAAAGGGTGCCCGA’, IUPACUnambiguousDNA())
>>> mutable_seq.remove("T")
>>> mutable_seq
MutableSeq(’GCCACGTAATGGGCCGCTGAAAGGGTGCCCGA’, IUPACUnambiguousDNA())
>>> mutable_seq.reverse()
>>> mutable_seq
MutableSeq(’AGCCCGTGGGAAAGTCGCCGGGTAATGCACCG’, IUPACUnambiguousDNA())
 
当然你也可以很容易的将上述生成的可变序列转换成只读模式：
>>> new_seq = mutable_seq.toseq()
>>> new_seq
Seq(’AGCCCGTGGGAAAGTCGCCGGGTAATGCACCG’, IUPACUnambiguousDNA())
 
 13.未知序列对象
 
当你有一段未知序列时候，你仅仅知道他的长度的时候，这个方法就用到了。
>>> from Bio.Seq import UnknownSeq
>>> unk = UnknownSeq(20)
>>> unk
UnknownSeq(20, alphabet = Alphabet(), character = ’?’)  #你也可以把？换成N等字符
>>> print unk
????????????????????
>>> len(unk)
20
 
你也可以针对性的处理未知序列：
>>> unk_dna
UnknownSeq(20, alphabet = IUPACAmbiguousDNA(), character = ’N’)
>>> unk_dna.complement()
UnknownSeq(20, alphabet = IUPACAmbiguousDNA(), character = ’N’)
>>> unk_dna.reverse_complement()
UnknownSeq(20, alphabet = IUPACAmbiguousDNA(), character = ’N’)
>>> unk_dna.transcribe()
UnknownSeq(20, alphabet = IUPACAmbiguousRNA(), character = ’N’)
>>> unk_protein = unk_dna.translate()
>>> unk_protein
UnknownSeq(6, alphabet = ProteinAlphabet(), character = ’X’)
>>> print unk_protein
XXXXXX
>>> len(unk_protein)
6

14.直接处理字符串
 
对于那些不想使用上述方法的读者，也可以直接对字符串进行处理的：
>>> 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*’
但是并不鼓励你用这种方法。
详见百度空间：http://hi.baidu.com/ilovelittree/item/750a4502e606ae03a1312de1