序列对象

生物序列可以说是生物信息学中的核心对象,在本节中,我们将介绍 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 字符串有非常有用的 upperlower 方法,用于更改大小写。例如,

>>> 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 片段,它编码一个短肽

\[\begin{split}\begin{gathered} \text{DNA 编码链(又称克里克链,链 } +1 \text{)} \\ \text{5'} \qquad \texttt{ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG} \qquad \text{3'} \\ \texttt{|||||||||||||||||||||||||||||||||||||||} \\ \text{3'} \qquad \texttt{TACCGGTAACATTACCCGGCGACTTTCCCACGGGCTATC} \qquad \text{5'} \\ \text{DNA 模板链(又称沃森链,链 } -1 \text{)} \end{gathered}\end{split}\]

此 DNA 序列的转录产生以下 RNA 序列

\[\begin{split}\begin{gathered} \text{5'} \qquad \texttt{AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG} \qquad \text{3'} \\ \text{单链信使 RNA} \end{gathered}\end{split}\]

实际的生物转录过程是从模板链开始,进行反向互补(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 对象的 transcribeback_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 序列,而没有明确指定序列内容。这种序列可以通过创建具有参数 NoneSeq 对象,然后是序列长度来表示。

>>> 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 对象转换为字符串 节)。

查找子序列

序列对象具有 findrfindindexrindex 方法,它们执行与普通字符串对象上的相应方法相同的函数。唯一的区别是子序列可以是字符串 (str)、bytesbytearraySeqMutableSeq 对象。

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

方法 rfindrindex 从序列的右侧开始搜索子序列。

>>> 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 方法还将普通字符串、bytesbytearraySeqMutableSeq 对象作为子序列;与上述示例一样,相同的子序列只报告一次。

直接使用字符串

为了结束本章,对于那些真的不想使用序列对象(或者更喜欢函数式编程风格而不是面向对象风格)的人,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 对象。