发那个太丢人

导航

 

转载
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

posted on 2018-05-26 14:33  发那个太丢人  阅读(598)  评论(0编辑  收藏  举报