序列对象
生物序列可以说是生物信息学中的核心对象,在本节中,我们将介绍 Biopython 处理序列的机制,即 Seq
对象。第 序列注释对象 将介绍相关的 SeqRecord
对象,它将序列信息与任何注释结合起来,在第 序列输入/输出 中用于序列输入/输出。
序列本质上是字母字符串,例如 AGTACACTGGT
,这看起来很自然,因为这是序列在生物文件格式中最常见的显示方式。
Seq
对象与标准 Python 字符串之间最重要的区别在于它们的 方法不同。虽然 Seq
对象支持与普通字符串相同的许多方法,但它的 translate()
方法通过执行生物学翻译而有所不同,并且还有一些其他生物学相关的 方法,例如 reverse_complement()
。
序列就像字符串
在大多数情况下,我们可以像处理普通 Python 字符串一样处理 Seq 对象,例如获取长度或迭代元素
>>> from Bio.Seq import Seq
>>> my_seq = Seq("GATCG")
>>> for index, letter in enumerate(my_seq):
... print("%i %s" % (index, letter))
...
0 G
1 A
2 T
3 C
4 G
>>> print(len(my_seq))
5
你可以像访问字符串的元素一样访问序列的元素(但请记住,Python 从零开始计数!)
>>> print(my_seq[0]) # first letter
G
>>> print(my_seq[2]) # third letter
T
>>> print(my_seq[-1]) # last letter
G
Seq
对象有一个 .count()
方法,就像字符串一样。请注意,这意味着就像 Python 字符串一样,这将给出非重叠计数
>>> from Bio.Seq import Seq
>>> "AAAA".count("AA")
2
>>> Seq("AAAA").count("AA")
2
对于某些生物学用途,你可能实际上想要一个重叠计数(即本例中的 \(3\))。当搜索单个字母时,这没有区别
>>> from Bio.Seq import Seq
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")
>>> len(my_seq)
32
>>> my_seq.count("G")
9
>>> 100 * (my_seq.count("G") + my_seq.count("C")) / len(my_seq)
46.875
虽然可以使用上面的代码片段来计算 GC%,但请注意 Bio.SeqUtils
模块已经内置了几个 GC 函数。例如
>>> from Bio.Seq import Seq
>>> from Bio.SeqUtils import gc_fraction
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")
>>> gc_fraction(my_seq)
0.46875
请注意,使用 Bio.SeqUtils.gc_fraction()
函数应该自动处理混合大小写的序列和模糊核苷酸 S,它代表 G 或 C。
另请注意,就像普通 Python 字符串一样,Seq
对象在某种程度上是“只读的”。如果你需要编辑序列,例如模拟点突变,请查看下面的第 MutableSeq 对象 节,其中介绍了 MutableSeq
对象。
切片序列
一个更复杂的例子,让我们获取序列的一部分
>>> from Bio.Seq import Seq
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")
>>> my_seq[4:12]
Seq('GATGGGCC')
请注意,“Seq” 对象遵循 Python 字符串的常规索引约定,序列的第一个元素编号为 0。当你进行切片时,第一个项目会被包含(即在本例中的 4),最后一个项目会被排除(在本例中的 12)。
同样像 Python 字符串一样,你可以使用起始、结束和步长(步长,默认为 1)来执行切片。例如,我们可以获取此 DNA 序列的第一、第二和第三密码子位置
>>> my_seq[0::3]
Seq('GCTGTAGTAAG')
>>> my_seq[1::3]
Seq('AGGCATGCATC')
>>> my_seq[2::3]
Seq('TAGCTAAGAC')
你可能在使用 Python 字符串时见过的另一个步长技巧是使用 -1 步长来反转字符串。你也可以用 Seq
对象来实现这一点
>>> my_seq[::-1]
Seq('CGCTAAAAGCTAGGATATATCCGGGTAGCTAG')
将 Seq 对象转换为字符串
如果你真的只需要一个普通字符串,例如写入文件或插入数据库,那么这很容易实现
>>> str(my_seq)
'GATCGATGGGCCTATATAGGATCGAAAATCGC'
由于对 Seq
对象调用 str()
会将整个序列作为字符串返回,因此你通常不需要显式地执行此转换。Python 在 print 函数中会自动执行此操作
>>> print(my_seq)
GATCGATGGGCCTATATAGGATCGAAAATCGC
你也可以在使用 Python 字符串格式化或插值运算符(%
)时,使用 Seq
对象直接使用 %s
占位符
>>> fasta_format_string = ">Name\n%s\n" % my_seq
>>> print(fasta_format_string)
>Name
GATCGATGGGCCTATATAGGATCGAAAATCGC
这行代码构建了一个简单的 FASTA 格式记录(无需担心换行)。第 format 方法 描述了一种从 SeqRecord
对象获取 FASTA 格式字符串的巧妙方法,而更普遍的关于读取和写入 FASTA 格式序列文件的话题将在第 序列输入/输出 中介绍。
连接或添加序列
两个 Seq
对象可以通过加法连接在一起
>>> from Bio.Seq import Seq
>>> seq1 = Seq("ACGT")
>>> seq2 = Seq("AACCGG")
>>> seq1 + seq2
Seq('ACGTAACCGG')
Biopython 不会检查序列内容,也不会在例如你连接蛋白质序列和 DNA 序列(这很可能是错误的)时引发异常
>>> from Bio.Seq import Seq
>>> protein_seq = Seq("EVRNAK")
>>> dna_seq = Seq("ACGT")
>>> protein_seq + dna_seq
Seq('EVRNAKACGT')
你可能经常需要将许多序列加在一起,这可以使用以下循环来实现
>>> from Bio.Seq import Seq
>>> list_of_seqs = [Seq("ACGT"), Seq("AACC"), Seq("GGTT")]
>>> concatenated = Seq("")
>>> for s in list_of_seqs:
... concatenated += s
...
>>> concatenated
Seq('ACGTAACCGGTT')
与 Python 字符串一样,Biopython Seq
也具有 .join
方法
>>> from Bio.Seq import Seq
>>> contigs = [Seq("ATG"), Seq("ATCCCG"), Seq("TTGCA")]
>>> spacer = Seq("N" * 10)
>>> spacer.join(contigs)
Seq('ATGNNNNNNNNNNATCCCGNNNNNNNNNNTTGCA')
更改大小写
Python 字符串有非常有用的 upper
和 lower
方法,用于更改大小写。例如,
>>> from Bio.Seq import Seq
>>> dna_seq = Seq("acgtACGT")
>>> dna_seq
Seq('acgtACGT')
>>> dna_seq.upper()
Seq('ACGTACGT')
>>> dna_seq.lower()
Seq('acgtacgt')
这些对于执行不区分大小写的匹配很有用
>>> "GTAC" in dna_seq
False
>>> "GTAC" in dna_seq.upper()
True
核苷酸序列和(反向)互补序列
对于核苷酸序列,你可以使用其内置方法轻松地获得 Seq
对象的互补序列或反向互补序列
>>> from Bio.Seq import Seq
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")
>>> my_seq
Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC')
>>> my_seq.complement()
Seq('CTAGCTACCCGGATATATCCTAGCTTTTAGCG')
>>> my_seq.reverse_complement()
Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC')
如前所述,反转 Seq
对象(或 Python 字符串)的简单方法是使用 -1 步长对其进行切片
>>> my_seq[::-1]
Seq('CGCTAAAAGCTAGGATATATCCGGGTAGCTAG')
如果你不小心尝试执行一些奇怪的操作,例如获取蛋白质序列的(反向)互补序列,结果将没有生物学意义
>>> from Bio.Seq import Seq
>>> protein_seq = Seq("EVRNAK")
>>> protein_seq.complement()
Seq('EBYNTM')
这里字母“E”不是核苷酸的有效 IUPAC 模糊代码,因此没有进行互补。但是,“V”表示“A”、“C”或“G”,其互补序列为“B”,依此类推。
第 将一组序列文件转换为它们的逆向互补序列 中的示例结合了 Seq
对象的逆向互补方法和 Bio.SeqIO
用于序列输入/输出。
转录
在讨论转录之前,我想澄清一下链问题。考虑以下(虚构的)双链 DNA 片段,它编码一个短肽
此 DNA 序列的转录产生以下 RNA 序列
实际的生物转录过程是从模板链开始,进行反向互补(TCAG \(→\) CUGA)以得到 mRNA。但是,在 Biopython 和生物信息学中,我们通常直接使用编码链,因为这意味着我们可以通过简单地将 T \(→\) U 来获得 mRNA 序列。
现在让我们真正开始在 Biopython 中执行转录。首先,让我们为编码链和模板 DNA 链创建 Seq
对象
>>> from Bio.Seq import Seq
>>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
>>> coding_dna
Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG')
>>> template_dna = coding_dna.reverse_complement()
>>> template_dna
Seq('CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT')
这些应该与上面的图一致 - 请记住,按照惯例,核苷酸序列通常从 5’ 到 3’ 方向读取,而在图中,模板链是反向显示的。
现在让我们使用 Seq
对象的内置 transcribe
方法将编码链转录成相应的 mRNA
>>> coding_dna
Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG')
>>> messenger_rna = coding_dna.transcribe()
>>> messenger_rna
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG')
正如你所见,这只是将 T 替换为 U。
如果你确实想要从模板链开始进行真正的生物转录,那么这将是一个两步过程
>>> template_dna.reverse_complement().transcribe()
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG')
Seq
对象还包含一个反向转录方法,用于从 mRNA 到 DNA 的编码链。同样,这是一个简单的 U \(→\) T 替换
>>> from Bio.Seq import Seq
>>> messenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG")
>>> messenger_rna
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG')
>>> messenger_rna.back_transcribe()
Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG')
注意:Seq
对象的 transcribe
和 back_transcribe
方法是在 Biopython 1.49 中添加的。对于较旧的版本,你需要使用 Bio.Seq
模块的函数,请查看第 直接使用字符串 节。
翻译
继续使用上面转录部分中讨论的相同示例,现在让我们将此 mRNA 翻译成相应的蛋白质序列 - 再次利用 Seq
对象的生物学方法之一
>>> from Bio.Seq import Seq
>>> messenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG")
>>> messenger_rna
Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG')
>>> messenger_rna.translate()
Seq('MAIVMGR*KGAR*')
您也可以直接从编码链 DNA 序列进行翻译。
>>> from Bio.Seq import Seq
>>> coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
>>> coding_dna
Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG')
>>> coding_dna.translate()
Seq('MAIVMGR*KGAR*')
您应该注意到上述蛋白质序列中,除了终止符之外,还有一个内部终止符。这是一个刻意选择的例子,因为它提供了一个解释一些可选参数的理由,包括不同的翻译表(遗传密码)。
Biopython 中提供的翻译表基于 NCBI 的那些(参见本教程的下一节)。默认情况下,翻译将使用标准遗传密码(NCBI 表格 ID 1)。假设我们正在处理一个线粒体序列。我们需要告诉翻译函数使用相关的遗传密码。
>>> coding_dna.translate(table="Vertebrate Mitochondrial")
Seq('MAIVMGRWKGAR*')
您也可以使用 NCBI 表格编号来指定表格,它更短,并且通常包含在 GenBank 文件的特征注释中。
>>> coding_dna.translate(table=2)
Seq('MAIVMGRWKGAR*')
现在,您可能想要将核苷酸翻译到第一个同框终止密码子,然后停止(就像在自然界中发生的那样)。
>>> coding_dna.translate()
Seq('MAIVMGR*KGAR*')
>>> coding_dna.translate(to_stop=True)
Seq('MAIVMGR')
>>> coding_dna.translate(table=2)
Seq('MAIVMGRWKGAR*')
>>> coding_dna.translate(table=2, to_stop=True)
Seq('MAIVMGRWKGAR')
请注意,当您使用 to_stop
参数时,终止密码子本身不会被翻译 - 并且终止符号不会包含在您蛋白质序列的末尾。
如果您不喜欢默认的星号,您甚至可以指定终止符号。
>>> coding_dna.translate(table=2, stop_symbol="@")
Seq('MAIVMGRWKGAR@')
现在,假设您有一个完整的编码序列 CDS,也就是说一个核苷酸序列(例如 mRNA - 经过任何剪接后),它是一个完整的密码子数(即长度是 3 的倍数),以起始密码子开始,以终止密码子结束,并且没有内部同框终止密码子。一般来说,给定一个完整的 CDS,默认的翻译方法会按您的要求进行(可能使用 to_stop
选项)。但是,如果您的序列使用了一个非标准的起始密码子呢?这在细菌中经常发生 - 例如,大肠杆菌 K12 中的 yaaX 基因。
>>> from Bio.Seq import Seq
>>> gene = Seq(
... "GTGAAAAAGATGCAATCTATCGTACTCGCACTTTCCCTGGTTCTGGTCGCTCCCATGGCA"
... "GCACAGGCTGCGGAAATTACGTTAGTCCCGTCAGTAAAATTACAGATAGGCGATCGTGAT"
... "AATCGTGGCTATTACTGGGATGGAGGTCACTGGCGCGACCACGGCTGGTGGAAACAACAT"
... "TATGAATGGCGAGGCAATCGCTGGCACCTACACGGACCGCCGCCACCGCCGCGCCACCAT"
... "AAGAAAGCTCCTCATGATCATCACGGCGGTCATGGTCCAGGCAAACATCACCGCTAA"
... )
>>> gene.translate(table="Bacterial")
Seq('VKKMQSIVLALSLVLVAPMAAQAAEITLVPSVKLQIGDRDNRGYYWDGGHWRDH...HR*',
ProteinAlpabet())
>>> gene.translate(table="Bacterial", to_stop=True)
Seq('VKKMQSIVLALSLVLVAPMAAQAAEITLVPSVKLQIGDRDNRGYYWDGGHWRDH...HHR')
在细菌遗传密码中,GTG
是一个有效的起始密码子,虽然它通常编码缬氨酸,但如果用作起始密码子,它应该被翻译成蛋氨酸。如果您告诉 Biopython 您的序列是一个完整的 CDS,就会发生这种情况。
>>> gene.translate(table="Bacterial", cds=True)
Seq('MKKMQSIVLALSLVLVAPMAAQAAEITLVPSVKLQIGDRDNRGYYWDGGHWRDH...HHR')
除了告诉 Biopython 将替代起始密码子翻译成蛋氨酸之外,使用此选项还可以确保您的序列确实是有效的 CDS(如果不是,您将得到一个异常)。
第 翻译 CDS 条目的 FASTA 文件 节中的示例将 Seq
对象的翻译方法与 Bio.SeqIO
结合起来用于序列输入/输出。
翻译表
在前面的章节中,我们讨论了 Seq
对象的翻译方法(并提到了 Bio.Seq
模块中的等效函数 - 参见第 直接使用字符串 节)。在内部,它们使用从 ftp://ftp.ncbi.nlm.nih.gov/entrez/misc/data/gc.prt 中的 NCBI 信息派生的密码子表对象,也在 https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi 上以更易读的布局显示。
和以前一样,我们只关注两个选择:标准翻译表和脊椎动物线粒体 DNA 的翻译表。
>>> from Bio.Data import CodonTable
>>> standard_table = CodonTable.unambiguous_dna_by_name["Standard"]
>>> mito_table = CodonTable.unambiguous_dna_by_name["Vertebrate Mitochondrial"]
或者,这些表格分别用 ID 号 1 和 2 标记。
>>> from Bio.Data import CodonTable
>>> standard_table = CodonTable.unambiguous_dna_by_id[1]
>>> mito_table = CodonTable.unambiguous_dna_by_id[2]
您可以通过打印它们来直观地比较实际的表格。
>>> print(standard_table)
Table 1 Standard, SGC0
| T | C | A | G |
--+---------+---------+---------+---------+--
T | TTT F | TCT S | TAT Y | TGT C | T
T | TTC F | TCC S | TAC Y | TGC C | C
T | TTA L | TCA S | TAA Stop| TGA Stop| A
T | TTG L(s)| TCG S | TAG Stop| TGG W | G
--+---------+---------+---------+---------+--
C | CTT L | CCT P | CAT H | CGT R | T
C | CTC L | CCC P | CAC H | CGC R | C
C | CTA L | CCA P | CAA Q | CGA R | A
C | CTG L(s)| CCG P | CAG Q | CGG R | G
--+---------+---------+---------+---------+--
A | ATT I | ACT T | AAT N | AGT S | T
A | ATC I | ACC T | AAC N | AGC S | C
A | ATA I | ACA T | AAA K | AGA R | A
A | ATG M(s)| ACG T | AAG K | AGG R | G
--+---------+---------+---------+---------+--
G | GTT V | GCT A | GAT D | GGT G | T
G | GTC V | GCC A | GAC D | GGC G | C
G | GTA V | GCA A | GAA E | GGA G | A
G | GTG V | GCG A | GAG E | GGG G | G
--+---------+---------+---------+---------+--
以及
>>> print(mito_table)
Table 2 Vertebrate Mitochondrial, SGC1
| T | C | A | G |
--+---------+---------+---------+---------+--
T | TTT F | TCT S | TAT Y | TGT C | T
T | TTC F | TCC S | TAC Y | TGC C | C
T | TTA L | TCA S | TAA Stop| TGA W | A
T | TTG L | TCG S | TAG Stop| TGG W | G
--+---------+---------+---------+---------+--
C | CTT L | CCT P | CAT H | CGT R | T
C | CTC L | CCC P | CAC H | CGC R | C
C | CTA L | CCA P | CAA Q | CGA R | A
C | CTG L | CCG P | CAG Q | CGG R | G
--+---------+---------+---------+---------+--
A | ATT I(s)| ACT T | AAT N | AGT S | T
A | ATC I(s)| ACC T | AAC N | AGC S | C
A | ATA M(s)| ACA T | AAA K | AGA Stop| A
A | ATG M(s)| ACG T | AAG K | AGG Stop| G
--+---------+---------+---------+---------+--
G | GTT V | GCT A | GAT D | GGT G | T
G | GTC V | GCC A | GAC D | GGC G | C
G | GTA V | GCA A | GAA E | GGA G | A
G | GTG V(s)| GCG A | GAG E | GGG G | G
--+---------+---------+---------+---------+--
您可能会发现以下属性很有用 - 例如,如果您尝试进行自己的基因查找。
>>> mito_table.stop_codons
['TAA', 'TAG', 'AGA', 'AGG']
>>> mito_table.start_codons
['ATT', 'ATC', 'ATA', 'ATG', 'GTG']
>>> mito_table.forward_table["ACG"]
'T'
比较 Seq 对象
序列比较实际上是一个非常复杂的话题,没有简单的方法来判断两个序列是否相等。基本问题是序列中字母的含义取决于上下文 - 字母“A”可能是 DNA、RNA 或蛋白质序列的一部分。Biopython 可以跟踪分子类型,因此比较两个 Seq
对象可能意味着也要考虑这一点。
DNA 片段“ACG”和 RNA 片段“ACG”应该相等吗?肽“ACG”呢?或者 Python 字符串“ACG”?在日常使用中,您的序列通常都是相同类型的(全部是 DNA、全部是 RNA,或者全部是蛋白质)。好吧,从 Biopython 1.65 开始,序列比较只查看序列,并像 Python 字符串一样进行比较。
>>> from Bio.Seq import Seq
>>> seq1 = Seq("ACGT")
>>> "ACGT" == seq1
True
>>> seq1 == "ACGT"
True
作为扩展,使用序列对象作为 Python 字典中的键等同于使用序列作为键的普通字符串。另请参见第 将 Seq 对象转换为字符串 节。
具有未知序列内容的序列
在某些情况下,可能知道序列的长度,但不知道构成它的实际字母。例如,GenBank 和 EMBL 文件可能只通过其配置信息来表示基因组 DNA 序列,而没有明确指定序列内容。这种序列可以通过创建具有参数 None
的 Seq
对象,然后是序列长度来表示。
>>> from Bio.Seq import Seq
>>> unknown_seq = Seq(None, 10)
因此创建的 Seq
对象具有明确定义的长度。但是,任何尝试访问序列内容的尝试都会引发 UndefinedSequenceError
。
>>> unknown_seq
Seq(None, length=10)
>>> len(unknown_seq)
10
>>> print(unknown_seq)
Traceback (most recent call last):
...
Bio.Seq.UndefinedSequenceError: Sequence content is undefined
>>>
具有部分定义序列内容的序列
有时,序列内容只对序列的部分进行定义,而在其他地方未定义。例如,MAF(多重比对格式)文件的以下摘录显示了人类、黑猩猩、猕猴、小鼠、大鼠、狗和负鼠基因组序列的比对。
s hg38.chr7 117512683 36 + 159345973 TTGAAAACCTGAATGTGAGAGTCAGTCAAGGATAGT
s panTro4.chr7 119000876 36 + 161824586 TTGAAAACCTGAATGTGAGAGTCACTCAAGGATAGT
s rheMac3.chr3 156330991 36 + 198365852 CTGAAATCCTGAATGTGAGAGTCAATCAAGGATGGT
s mm10.chr6 18207101 36 + 149736546 CTGAAAACCTAAGTAGGAGAATCAACTAAGGATAAT
s rn5.chr4 42326848 36 + 248343840 CTGAAAACCTAAGTAGGAGAGACAGTTAAAGATAAT
s canFam3.chr14 56325207 36 + 60966679 TTGAAAAACTGATTATTAGAGTCAATTAAGGATAGT
s monDom5.chr8 173163865 36 + 312544902 TTAAGAAACTGGAAATGAGGGTTGAATGACAAACTT
在每一行中,第一个数字表示比对序列在染色体上的起始位置(以零为基准的坐标),后面跟着比对序列的大小、链、完整染色体的大小以及比对序列。
可以使用字典作为 data
参数来创建表示这种部分定义序列的 Seq
对象,其中键是已知序列段的起始坐标,而值是相应的序列内容。例如,对于第一个序列,我们将使用
>>> from Bio.Seq import Seq
>>> seq = Seq({117512683: "TTGAAAACCTGAATGTGAGAGTCAGTCAAGGATAGT"}, length=159345973)
从部分定义序列中提取子序列可能会返回一个完全定义的序列、一个未定义的序列或一个部分定义的序列,具体取决于坐标。
>>> seq[1000:1020]
Seq(None, length=20)
>>> seq[117512690:117512700]
Seq('CCTGAATGTG')
>>> seq[117512670:117512690]
Seq({13: 'TTGAAAA'}, length=20)
>>> seq[117512700:]
Seq({0: 'AGAGTCAGTCAAGGATAGT'}, length=41833273)
如果至少有一个序列是部分或完全未定义的,则也可以通过追加序列来创建部分定义的序列。
>>> seq = Seq("ACGT")
>>> undefined_seq = Seq(None, length=10)
>>> seq + undefined_seq + seq
Seq({0: 'ACGT', 14: 'ACGT'}, length=18)
MutableSeq 对象
就像普通的 Python 字符串一样,Seq
对象是“只读”的,或者用 Python 术语来说,是不可变的。除了想要 Seq
对象的行为像一个字符串之外,这在许多生物学应用中也是一个有用的默认值,因为您希望确保没有更改您的序列数据。
>>> from Bio.Seq import Seq
>>> my_seq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")
观察一下如果您尝试编辑序列会发生什么。
>>> my_seq[5] = "G"
Traceback (most recent call last):
...
TypeError: 'Seq' object does not support item assignment
但是,您可以将其转换为可变序列(一个 MutableSeq
对象),然后就可以对其执行您想要的任何操作。
>>> from Bio.Seq import MutableSeq
>>> mutable_seq = MutableSeq(my_seq)
>>> mutable_seq
MutableSeq('GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA')
或者,您可以直接从字符串创建一个 MutableSeq
对象。
>>> from Bio.Seq import MutableSeq
>>> mutable_seq = MutableSeq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")
无论哪种方式,您都将获得一个可以更改的序列对象。
>>> mutable_seq
MutableSeq('GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA')
>>> mutable_seq[5] = "C"
>>> mutable_seq
MutableSeq('GCCATCGTAATGGGCCGCTGAAAGGGTGCCCGA')
>>> mutable_seq.remove("T")
>>> mutable_seq
MutableSeq('GCCACGTAATGGGCCGCTGAAAGGGTGCCCGA')
>>> mutable_seq.reverse()
>>> mutable_seq
MutableSeq('AGCCCGTGGGAAAGTCGCCGGGTAATGCACCG')
请注意,MutableSeq
对象的 reverse()
方法,就像 Python 列表的 reverse()
方法一样,会就地反转序列。
Python 中可变对象和不可变对象之间的一个重要技术差异意味着您不能将 MutableSeq
对象用作字典键,但您可以以这种方式使用 Python 字符串或 Seq
对象。
完成对 MutableSeq
对象的编辑后,如果需要,可以轻松地将其转换回只读的 Seq
对象。
>>> from Bio.Seq import Seq
>>> new_seq = Seq(mutable_seq)
>>> new_seq
Seq('AGCCCGTGGGAAAGTCGCCGGGTAATGCACCG')
您也可以从 MutableSeq
对象获取字符串,就像从 Seq
对象获取字符串一样(参见第 将 Seq 对象转换为字符串 节)。
查找子序列
序列对象具有 find
、rfind
、index
和 rindex
方法,它们执行与普通字符串对象上的相应方法相同的函数。唯一的区别是子序列可以是字符串 (str
)、bytes
、bytearray
、Seq
或 MutableSeq
对象。
>>> from Bio.Seq import Seq, MutableSeq
>>> seq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")
>>> seq.index("ATGGGCCGC")
9
>>> seq.index(b"ATGGGCCGC")
9
>>> seq.index(bytearray(b"ATGGGCCGC"))
9
>>> seq.index(Seq("ATGGGCCGC"))
9
>>> seq.index(MutableSeq("ATGGGCCGC"))
9
如果找不到子序列,则会引发 ValueError
。
>>> seq.index("ACTG")
Traceback (most recent call last):
...
ValueError: ...
而 find
方法如果找不到子序列,则返回 -1。
>>> seq.find("ACTG")
-1
方法 rfind
和 rindex
从序列的右侧开始搜索子序列。
>>> seq.find("CC")
1
>>> seq.rfind("CC")
29
使用 search
方法同时搜索多个子序列。此方法返回一个迭代器。
>>> for index, sub in seq.search(["CC", "GGG", "CC"]):
... print(index, sub)
...
1 CC
11 GGG
14 CC
23 GGG
28 CC
29 CC
search
方法还将普通字符串、bytes
、bytearray
、Seq
和 MutableSeq
对象作为子序列;与上述示例一样,相同的子序列只报告一次。
直接使用字符串
为了结束本章,对于那些真的不想使用序列对象(或者更喜欢函数式编程风格而不是面向对象风格)的人,Bio.Seq
中的模块级函数将接受普通的 Python 字符串、Seq
对象或 MutableSeq
对象。
>>> 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*'
但是,鼓励您默认使用 Seq
对象。