序列注释对象
第 序列对象章介绍了序列类。直接位于 Seq
类之上的是序列记录或 SeqRecord
类,定义在 Bio.SeqRecord
模块中。此类允许将更高级别的功能(如标识符和特征(作为 SeqFeature
对象))与序列关联,并且在整个序列输入/输出接口 Bio.SeqIO
中使用,该接口在第 序列输入/输出章中进行了详细介绍。
如果你只打算处理简单的 FASTA 文件等数据,你现在就可以跳过本章。另一方面,如果你要使用注释丰富的序列数据,例如来自 GenBank 或 EMBL 文件的数据,这些信息非常重要。
虽然本章应该涵盖本章中与 SeqRecord
和 SeqFeature
对象相关的多数操作,但你可能还想阅读 SeqRecord
wiki 页面(https://biopython.pythonlang.cn/wiki/SeqRecord)和内置文档(Bio.SeqRecord
和 Bio.SeqFeature
)。
>>> from Bio.SeqRecord import SeqRecord
>>> help(SeqRecord)
SeqRecord 对象
SeqRecord(序列记录)类定义在 Bio.SeqRecord
模块中。此类允许将更高级别的功能(如标识符和特征)与序列关联(见第 序列对象章),并且是 Bio.SeqIO
序列输入/输出接口的基本数据类型(见第 序列输入/输出章)。
SeqRecord 类本身非常简单,并以属性的形式提供以下信息:
- .seq
序列本身,通常是
Seq
对象。- .id
用于标识序列的主要 ID - 字符串。在大多数情况下,这类似于登录号。
- .name
序列的“常用”名称/ID - 字符串。在某些情况下,这将与登录号相同,但也可能是克隆名称。我认为这类似于 GenBank 记录中的 LOCUS ID。
- .description
序列的人类可读描述或表达名称 - 字符串。
- .letter_annotations
使用(受限的)字典存储每个字母的注释,这些注释包含有关序列中字母的附加信息。键是信息的名称,信息包含在值中,值是 Python 序列(即列表、元组或字符串),其长度与序列本身相同。这通常用于质量评分(例如,第 FASTQ 文件的简单质量过滤节)或二级结构信息(例如,来自 Stockholm/PFAM 比对文件)。
- .annotations
关于序列的附加信息的字典。键是信息的名称,信息包含在值中。这允许向序列添加更多“非结构化”信息。
- .features
SeqFeature 对象的列表,其中包含有关序列上的特征的更多结构化信息(例如,基因在基因组中的位置或蛋白质序列中的域)。序列特征的结构将在本章的第 特征、位置和位置对象节中进行描述。
- .dbxrefs
数据库交叉引用的字符串列表。
创建 SeqRecord
使用 SeqRecord 对象并不复杂,因为所有信息都以类的属性形式呈现。通常,你不会“手动”创建 SeqRecord,而是使用 Bio.SeqIO 来为你读取序列文件(见第 序列输入/输出章和下面的示例)。但是,创建 SeqRecord 非常简单。
从头开始创建 SeqRecord 对象
要创建 SeqRecord,至少需要一个 Seq 对象
>>> from Bio.Seq import Seq
>>> simple_seq = Seq("GATC")
>>> from Bio.SeqRecord import SeqRecord
>>> simple_seq_r = SeqRecord(simple_seq)
此外,你还可以将 ID、名称和描述传递给初始化函数,但如果没有,它们将被设置为字符串,表明它们是未知的,并且可以在之后修改
>>> simple_seq_r.id
'<unknown id>'
>>> simple_seq_r.id = "AC12345"
>>> simple_seq_r.description = "Made up sequence I wish I could write a paper about"
>>> print(simple_seq_r.description)
Made up sequence I wish I could write a paper about
>>> simple_seq_r.seq
Seq('GATC')
如果要将 SeqRecord 输出到文件,包含标识符非常重要。你通常会在创建对象时包含它
>>> from Bio.Seq import Seq
>>> simple_seq = Seq("GATC")
>>> from Bio.SeqRecord import SeqRecord
>>> simple_seq_r = SeqRecord(simple_seq, id="AC12345")
如上所述,SeqRecord 具有一个名为 annotations 的字典属性。它用于任何不适合其他更特定属性的杂项注释。添加注释很简单,只需直接处理注释字典即可
>>> simple_seq_r.annotations["evidence"] = "None. I just made it up."
>>> print(simple_seq_r.annotations)
{'evidence': 'None. I just made it up.'}
>>> print(simple_seq_r.annotations["evidence"])
None. I just made it up.
处理每个字母的注释类似,letter_annotations 是一个类似字典的属性,它允许你为任何 Python 序列(即字符串、列表或元组)赋值,其长度与序列相同
>>> simple_seq_r.letter_annotations["phred_quality"] = [40, 40, 38, 30]
>>> print(simple_seq_r.letter_annotations)
{'phred_quality': [40, 40, 38, 30]}
>>> print(simple_seq_r.letter_annotations["phred_quality"])
[40, 40, 38, 30]
dbxrefs 和 features 属性只是 Python 列表,应分别用于存储字符串和 SeqFeature 对象(本章稍后讨论)。
从 FASTA 文件创建 SeqRecord 对象
此示例使用一个相当大的 FASTA 文件,其中包含来自 NCBI 下载的鼠疫耶尔森氏菌生物变种微鼠型菌株 91001 质粒 pPCP1 的完整序列。此文件包含在 Biopython 单元测试中,位于 GenBank 文件夹下,或在线 NC_005816.fna(来自我们的网站)。
该文件以以下方式开头 - 并且你可以检查是否只有一个记录存在(即,只有一行以大于号开头)
>gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus ... pPCP1, complete sequence
TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGGGGGTAATCTGCTCTCC
...
在第 快速入门 - 你可以用 Biopython 做什么?章中,你将看到函数 Bio.SeqIO.parse(...) 用于循环遍历文件中所有的记录,这些记录作为 SeqRecord 对象。Bio.SeqIO 模块具有一个姐妹函数,用于处理只包含一个记录的文件,我们将在这里使用它(有关详细信息,请参阅第 序列输入/输出章)。
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.fna", "fasta")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='gi|45478711|ref|NC_005816.1|', name='gi|45478711|ref|NC_005816.1|', description='gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])
现在,让我们看一下此 SeqRecord 的关键属性 - 从 seq 属性开始,它提供了一个 Seq 对象
>>> record.seq
Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG')
接下来是标识符和描述
>>> record.id
'gi|45478711|ref|NC_005816.1|'
>>> record.name
'gi|45478711|ref|NC_005816.1|'
>>> record.description
'gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
如上所示,FASTA 记录标题行的第一个单词(删除大于号后)用于 id 和 name 属性。整个标题行(删除大于号后)用于记录描述。这是故意的,部分原因是出于向后兼容性的原因,但如果你有这样的 FASTA 文件,它也有道理
>Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1
TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGGGGGTAATCTGCTCTCC
...
请注意,读取 FASTA 文件时,其他注释属性都没有填充
>>> record.dbxrefs
[]
>>> record.annotations
{}
>>> record.letter_annotations
{}
>>> record.features
[]
在本例中,我们的示例 FASTA 文件来自 NCBI,他们有一套相当完善的格式化 FASTA 行的约定。这意味着可以解析此信息并提取 GI 编号和登录号等信息。但是,来自其他来源的 FASTA 文件各不相同,因此这在一般情况下是不可能的。
从 GenBank 文件创建 SeqRecord 对象
与之前的示例类似,我们将查看Yersinia pestis biovar Microtus菌株 91001 质粒 pPCP1 的完整序列,该序列最初从 NCBI 下载,但这次是作为 GenBank 文件。同样,该文件包含在 Biopython 单元测试中,位于 GenBank 文件夹下,或在我们的网站上在线 NC_005816.gb。
该文件包含单个记录(即只有一个 LOCUS 行)并以
LOCUS NC_005816 9609 bp DNA circular BCT 21-JUL-2008
DEFINITION Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete
sequence.
ACCESSION NC_005816
VERSION NC_005816.1 GI:45478711
PROJECT GenomeProject:10638
...
同样,我们将使用 Bio.SeqIO
读取此文件,代码与上面用于 FASTA 文件的代码几乎相同(有关详细信息,请参见第 序列输入/输出 章)
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])
>>> record.seq
Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG')
The name
来自 LOCUS 行,而 id
包括版本后缀。描述来自 DEFINITION 行
>>> record.id
'NC_005816.1'
>>> record.name
'NC_005816'
>>> record.description
'Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
GenBank 文件没有任何逐字母注释
>>> record.letter_annotations
{}
大多数注释信息都记录在 annotations
字典中,例如
>>> len(record.annotations)
13
>>> record.annotations["source"]
'Yersinia pestis biovar Microtus str. 91001'
The dbxrefs
列表从任何 PROJECT 或 DBLINK 行填充
>>> record.dbxrefs
['Project:58037']
最后,也许最有趣的是,特征表中的所有条目(例如基因或 CDS 特征)都作为 SeqFeature
对象记录在 features
列表中。
>>> len(record.features)
41
我们将在下一节 特征、位置和位置对象 中讨论 SeqFeature
对象。
特征、位置和位置对象
SeqFeature 对象
序列特征是描述序列的重要组成部分。一旦你超越了序列本身,你就需要某种方法来组织并轻松获取有关序列的更多“抽象”信息。虽然开发一个通用的序列特征类来涵盖所有内容可能是不可能的,但 Biopython SeqFeature
类试图封装尽可能多的有关序列的信息。该设计主要基于 GenBank/EMBL 特征表,因此,如果你理解它们的外观,你可能会更容易理解 Biopython 类别的结构。
每个 SeqFeature
对象的关键思想是描述父序列上的一个区域,通常是一个 SeqRecord
对象。该区域用一个位置对象描述,通常是两个位置之间的范围(见下文 位置和位置 部分)。
The SeqFeature
类具有许多属性,因此我们将首先列出它们及其通用特征,然后在本章的后面部分通过示例来展示如何在现实生活中应用它们。SeqFeature 的属性是
- .type
这是对特征类型进行文本描述(例如,这将是类似于“CDS”或“gene”的东西)。
- .location
The
SeqFeature
在你正在处理的序列上的位置,见下文 位置和位置 部分。TheSeqFeature
将其大部分功能委派给位置对象,并包括许多用于位置属性的快捷属性- .ref
是
.location.ref
的简写 - 位置引用的任何(不同的)参考序列。通常只是 None。- .ref_db
是
.location.ref_db
的简写 - 指定.ref
中的任何标识符所引用的数据库。通常只是 None。- .strand
是
.location.strand
的简写 - 特征位于序列上的链。对于双链核苷酸序列,这可能是 \(1\) 表示顶链,\(-1\) 表示底链,\(0\) 表示链很重要但未知,或者None
表示无关紧要。对于蛋白质或单链序列,这为 None。
- .qualifiers
这是一个 Python 字典,包含有关该特征的附加信息。键是某种简短的单字描述,说明值中包含的信息是什么,而值是实际信息。例如,限定符的常见键可能是“证据”,而值可能是“计算(非实验)”。这只是为了让查看该特征的人知道它尚未通过实验(即在湿实验室中)得到证实。请注意,其他值将是字符串列表(即使只有一个字符串)。这反映了 GenBank/EMBL 文件中的特征表。
- .sub_features
这曾经用于表示具有复杂位置的特征,例如 GenBank/EMBL 文件中的“连接”。随着
CompoundLocation
对象的引入,这已被弃用,现在应该忽略。
位置和位置
每个 SeqFeature
对象的关键思想是描述父序列上的一个区域,为此我们使用一个位置对象,通常描述两个位置之间的范围。有两个尝试澄清我们使用的术语
- 位置
这指的是序列上的单个位置,它可能是模糊的也可能不是。例如,5、20、
<100
和>200
都是位置。- 位置
位置是序列的区域,以一些位置为界。例如
5..20
(即 5 到 20)是一个位置。
我只是提到这一点,因为有时我会在这两者之间感到困惑。
SimpleLocation 对象
除非你处理真核基因,否则大多数 SeqFeature
位置都非常简单 - 你只需要起点和终点坐标以及一条链。这基本上是基本的 SimpleLocation
对象所做的全部工作。
在实践中,情况可能更复杂。首先,我们必须处理由几个区域组成的复合位置。其次,位置本身可能是模糊的(不精确)。
CompoundLocation 对象
Biopython 1.62 引入了 CompoundLocation
作为重构由多个区域组成的复杂位置表示方式的一部分。主要用途是处理 EMBL/GenBank 文件中的“连接”位置。
模糊位置
到目前为止,我们只使用了简单的位置。处理特征位置的一个复杂之处在于位置本身。在生物学中,很多时候事情并不完全确定(尽管我们湿实验室生物学家尽力使它们确定!)。例如,你可能会进行二核苷酸引物实验,并发现 mRNA 转录本的起始位置在两个位点中的一个。这是一个非常有用的信息,但复杂之处在于如何将它表示为一个位置。为了帮助我们处理这个问题,我们有模糊位置的概念。基本上有几种类型的模糊位置,因此我们有五个类来处理它们
- ExactPosition
顾名思义,此类表示一个位置,该位置被指定为沿着序列的精确位置。这被表示为一个数字,你可以通过查看对象的
position
属性来获取该位置。- BeforePosition
此类表示一个模糊位置,该位置发生在某个指定位点之前。在 GenBank/EMBL 符号中,这表示为类似于
<13
的东西,表示实际位置位于 13 以下的某个位置。要获取指定的边界,请查看对象的position
属性。- AfterPosition
与
BeforePosition
相反,此类表示一个位置,该位置发生在某个指定位点之后。这在 GenBank 中表示为>13
,与BeforePosition
类似,你可以通过查看对象的position
属性来获取边界数字。- WithinPosition
偶尔用于 GenBank/EMBL 位置,此类模拟一个位置,该位置发生在两个指定核苷酸之间。在 GenBank/EMBL 符号中,这将表示为
(1.5)
,以表示该位置位于 1 到 5 的范围内。- OneOfPosition
偶尔用于 GenBank/EMBL 位置,此类处理一个位置,其中存在多个可能的值,例如,如果起始密码子不清楚并且存在两个基因起始候选,则可以使用此类。或者,那可以明确地处理为两个相关的基因特征。
- UnknownPosition
此类处理位置未知的位置。这在 GenBank/EMBL 中未使用,但对应于 UniProt 中使用的“?”特征坐标。
以下是一个创建具有模糊端点的位置的示例
>>> from Bio import SeqFeature
>>> start_pos = SeqFeature.AfterPosition(5)
>>> end_pos = SeqFeature.BetweenPosition(9, left=8, right=9)
>>> my_location = SeqFeature.SimpleLocation(start_pos, end_pos)
请注意,Biopython 1.59 中一些模糊位置的细节发生了变化,特别是对于 BetweenPosition 和 WithinPosition,你现在必须明确指定哪个整数位置应用于切片等。对于起始位置,这通常是较低(左侧)的值,而对于结束位置,这通常是较高(右侧)的值。
如果你打印出 SimpleLocation
对象,你可以获得该信息的良好表示
>>> print(my_location)
[>5:(8^9)]
我们可以使用位置对象的 start 和 end 属性来访问模糊的起始位置和结束位置
>>> my_location.start
AfterPosition(5)
>>> print(my_location.start)
>5
>>> my_location.end
BetweenPosition(9, left=8, right=9)
>>> print(my_location.end)
(8^9)
如果你不想处理模糊位置,而只想获得数字,它们实际上是整数的子类,因此应该像整数一样工作
>>> int(my_location.start)
5
>>> int(my_location.end)
9
类似地,为了便于创建位置而不必担心模糊位置,你可以将数字传递给 FeaturePosition
构造函数,你将获得 ExactPosition
对象
>>> exact_location = SeqFeature.SimpleLocation(5, 9)
>>> print(exact_location)
[5:9]
>>> exact_location.start
ExactPosition(5)
>>> int(exact_location.start)
5
这就是在 Biopython 中处理模糊位置的大部分细节。它是经过设计的,因此处理模糊性并不比处理精确位置复杂多少,希望你发现这是真的!
位置测试
你可以将 Python 关键字 in
与 SeqFeature
或位置对象一起使用,以查看父坐标的碱基/残基是否在特征/位置内。
例如,假设您有一个感兴趣的 SNP,并且您想知道此 SNP 在哪些特征内,并且假设此 SNP 位于索引 4350(Python 计数!)。以下是一个简单的蛮力解决方案,我们只需在循环中逐个检查所有特征。
>>> from Bio import SeqIO
>>> my_snp = 4350
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> for feature in record.features:
... if my_snp in feature:
... print("%s %s" % (feature.type, feature.qualifiers.get("db_xref")))
...
source ['taxon:229193']
gene ['GeneID:2767712']
CDS ['GI:45478716', 'GeneID:2767712']
请注意,来自 GenBank 或 EMBL 文件的基因和 CDS 特征定义为连接,是外显子的并集 - 它们不覆盖任何内含子。
特征或位置描述的序列
一个 SeqFeature
或位置对象不直接包含序列,而是位置(参见第 位置和位置 节)描述了如何从父序列中获取它。例如,考虑一个具有位置 5:18
在反向链上的(短)基因序列,在 GenBank/EMBL 表示法中使用 1 为基的计数将是 complement(6..18)
,如下所示
>>> from Bio.Seq import Seq
>>> from Bio.SeqFeature import SeqFeature, SimpleLocation
>>> seq = Seq("ACCGAGACGGCAAAGGCTAGCATAGGTATGAGACTTCCTTCCTGCCAGTGCTGAGGAACTGGGAGCCTAC")
>>> feature = SeqFeature(SimpleLocation(5, 18, strand=-1), type="gene")
您可以获取父序列,对其进行切片以提取 5:18
,然后获取反向互补序列。特征位置的开始和结束是类似整数的,因此可以正常工作
>>> feature_seq = seq[feature.location.start : feature.location.end].reverse_complement()
>>> print(feature_seq)
AGCCTTTGCCGTC
这是一个简单的示例,因此这并不太糟糕 - 但是,一旦您必须处理复合特征(连接),这就会变得很乱。相反,SeqFeature
对象有一个 extract
方法来处理所有这些(并且从 Biopython 1.78 开始,可以通过提供一个引用序列的字典来处理跨剪接)
>>> feature_seq = feature.extract(seq)
>>> print(feature_seq)
AGCCTTTGCCGTC
一个 SeqFeature
或位置的长度与其描述的序列区域的长度相匹配。
>>> print(len(feature_seq))
13
>>> print(len(feature))
13
>>> print(len(feature.location))
13
对于 SimpleLocation
对象,长度只是开始位置和结束位置之间的差值。但是,对于 CompoundLocation
,长度是组成区域的总和。
比较
The SeqRecord
对象可能非常复杂,但这是一个简单的示例
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> record1 = SeqRecord(Seq("ACGT"), id="test")
>>> record2 = SeqRecord(Seq("ACGT"), id="test")
当您尝试比较这些“相同”记录时会发生什么?
>>> record1 == record2
也许令人惊讶的是,较旧版本的 Biopython 会为 SeqRecord
使用 Python 的默认对象比较,这意味着 record1 == record2
仅当这些变量指向内存中的同一个对象时才会返回 True
。在这个例子中,record1 == record2
会返回 False
!
>>> record1 == record2 # on old versions of Biopython!
False
从 Biopython 1.67 开始,SeqRecord
比较(如 record1 == record2
)会引发显式错误,以避免人们被此问题困扰
>>> record1 == record2
Traceback (most recent call last):
...
NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest.
相反,您应该检查您感兴趣的属性,例如标识符和序列
>>> record1.id == record2.id
True
>>> record1.seq == record2.seq
True
注意,比较复杂对象很快就会变得复杂(另见第 比较 Seq 对象 节)。
参考文献
另一个与序列相关的常见注释是对处理序列的期刊或其他已发表作品的引用。我们在 Biopython 中有一个相当简单的方法来表示引用 - 我们有一个 Bio.SeqFeature.Reference
类来存储有关引用的相关信息作为对象的属性。
这些属性包括您期望在引用中看到的属性,例如 journal
、title
和 authors
。此外,它还可以保存 medline_id
和 pubmed_id
以及关于引用的 comment
。所有这些都可以简单地作为对象的属性访问。
引用也具有 location
对象,以便它可以指定引用所指序列上的特定位置。例如,您可能有一篇关于位于 BAC 上的特定基因的期刊,并且希望指定它只引用此确切位置。The location
是一个可能模糊的位置,如第 位置和位置 节中所述。
任何引用对象都作为列表存储在 SeqRecord
对象的 annotations
字典中,键为“references”。这就是全部内容。引用应该易于处理,并且希望足够通用以涵盖许多使用案例。
格式化方法
The format()
方法 SeqRecord
类给出了一个字符串,其中包含使用 Bio.SeqIO
支持的输出文件格式(例如 FASTA)格式化的记录
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> record = SeqRecord(
... Seq(
... "MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD"
... "GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK"
... "NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM"
... "SSAC"
... ),
... id="gi|14150838|gb|AAK54648.1|AF376133_1",
... description="chalcone synthase [Cucumis sativus]",
... )
>>> print(record.format("fasta"))
这应该给出
>gi|14150838|gb|AAK54648.1|AF376133_1 chalcone synthase [Cucumis sativus]
MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD
GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK
NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM
SSAC
This format
方法接受一个必须的参数,一个由 Bio.SeqIO
作为输出格式支持的小写字符串(参见第 序列输入/输出 章)。但是,一些 Bio.SeqIO
可以写入的文件格式需要多条记录(通常是多序列比对格式的情况),因此无法通过此 format()
方法工作。另见第 将 SeqRecord 对象作为格式化的字符串获取 节。
切片 SeqRecord
您可以切片一个 SeqRecord
,以获得一个新的 SeqRecord
,该记录仅涵盖序列的一部分。这里重要的是,任何按字母的注释也会被切片,并且任何完全落在新序列内的特征都会被保留(其位置会调整)。
例如,使用之前使用的相同 GenBank 文件
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])
>>> len(record)
9609
>>> len(record.features)
41
对于此示例,我们将重点关注 pim
基因,YP_pPCP05
。如果您直接查看 GenBank 文件,您会发现此基因/CDS 的位置字符串为 4343..4780
,或在 Python 计数中为 4342:4780
。从查看该文件,您可以算出这些是文件中第 12 和第 13 个条目,因此在 Python 零为基的计数中,它们是 features
列表中的第 \(11\) 和 \(12\) 个条目
>>> print(record.features[20])
type: gene
location: [4342:4780](+)
qualifiers:
Key: db_xref, Value: ['GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
>>> print(record.features[21])
type: CDS
location: [4342:4780](+)
qualifiers:
Key: codon_start, Value: ['1']
Key: db_xref, Value: ['GI:45478716', 'GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
Key: note, Value: ['similar to many previously sequenced pesticin immunity protein entries of Yersinia pestis plasmid pPCP, e.g. gi| 16082683|,ref|NP_395230.1| (NC_003132) , gi|1200166|emb|CAA90861.1| (Z54145 ) , gi|1488655| emb|CAA63439.1| (X92856) , gi|2996219|gb|AAC62543.1| (AF053945) , and gi|5763814|emb|CAB531 67.1| (AL109969)']
Key: product, Value: ['pesticin immunity protein']
Key: protein_id, Value: ['NP_995571.1']
Key: transl_table, Value: ['11']
Key: translation, Value: ['MGGGMISKLFCLALIFLSSSGLAEKNTYTAKDILQNLELNTFGNSLSHGIYGKQTTFKQTEFTNIKSNTKKHIALINKDNSWMISLKILGIKRDEYTVCFEDFSLIRPPTYVAIHPLLIKKVKSGNFIVVKEIKKSIPGCTVYYH']
让我们从 4300 到 4800 切片此父记录(足以包含 pim
基因/CDS),看看我们获得了多少个特征
>>> sub_record = record[4300:4800]
>>> sub_record
SeqRecord(seq=Seq('ATAAATAGATTATTCCAAATAATTTATTTATGTAAGAACAGGATGGGAGGGGGA...TTA'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])
>>> len(sub_record)
500
>>> len(sub_record.features)
2
我们的子记录只有两个特征,即 YP_pPCP05
的基因和 CDS 条目
>>> print(sub_record.features[0])
type: gene
location: [42:480](+)
qualifiers:
Key: db_xref, Value: ['GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
>>> print(sub_record.features[1])
type: CDS
location: [42:480](+)
qualifiers:
Key: codon_start, Value: ['1']
Key: db_xref, Value: ['GI:45478716', 'GeneID:2767712']
Key: gene, Value: ['pim']
Key: locus_tag, Value: ['YP_pPCP05']
Key: note, Value: ['similar to many previously sequenced pesticin immunity protein entries of Yersinia pestis plasmid pPCP, e.g. gi| 16082683|,ref|NP_395230.1| (NC_003132) , gi|1200166|emb|CAA90861.1| (Z54145 ) , gi|1488655| emb|CAA63439.1| (X92856) , gi|2996219|gb|AAC62543.1| (AF053945) , and gi|5763814|emb|CAB531 67.1| (AL109969)']
Key: product, Value: ['pesticin immunity protein']
Key: protein_id, Value: ['NP_995571.1']
Key: transl_table, Value: ['11']
Key: translation, Value: ['MGGGMISKLFCLALIFLSSSGLAEKNTYTAKDILQNLELNTFGNSLSHGIYGKQTTFKQTEFTNIKSNTKKHIALINKDNSWMISLKILGIKRDEYTVCFEDFSLIRPPTYVAIHPLLIKKVKSGNFIVVKEIKKSIPGCTVYYH']
请注意,它们的位置已调整以反映新的父序列!
虽然 Biopython 已经对特征(以及任何按字母的注释)做了明智且直观的处理,但对于其他注释,我们无法知道它是否仍然适用于子序列。为了避免猜测,除了分子类型外,.annotations
和 .dbxrefs
将从子记录中省略,您需要根据需要显式传输任何相关信息。
>>> sub_record.annotations
{'molecule_type': 'DNA'}
>>> sub_record.dbxrefs
[]
您可能希望保留其他条目,例如生物体?注意复制整个 annotations 字典,因为在这种情况下,您的部分序列不再是环状 DNA - 它现在是线性序列
>>> sub_record.annotations["topology"] = "linear"
对于记录 id
、name
和 description
也可以这样说,但出于实用性考虑,这些将被保留
>>> sub_record.id
'NC_005816.1'
>>> sub_record.name
'NC_005816'
>>> sub_record.description
'Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
这很好地说明了这个问题,我们的新子记录不是质粒的完整序列,因此描述是错误的!让我们修复它,然后使用上面第 格式化方法 节中描述的 format
方法将子记录视为简化的 GenBank 文件
>>> sub_record.description = (
... "Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, partial"
... )
>>> print(sub_record.format("genbank")[:200] + "...")
LOCUS NC_005816 500 bp DNA linear UNK 01-JAN-1980
DEFINITION Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, partial.
ACCESSION NC_005816
VERSION NC_0058...
参见第 修剪掉引物序列 节和第 修剪掉接头序列 节,了解一些 FASTQ 示例,其中按字母的注释(读数质量分数)也会被切片。
添加 SeqRecord 对象
您可以将 SeqRecord
对象加在一起,得到一个新的 SeqRecord
。这里重要的是,任何常见的按字母的注释也会被添加,所有特征都会被保留(其位置会调整),并且任何其他常见的注释也会被保留(例如 id、name 和 description)。
对于带有按字母注释的示例,我们将使用 FASTQ 文件中的第一条记录。第 序列输入/输出 章将解释 SeqIO
函数
>>> from Bio import SeqIO
>>> record = next(SeqIO.parse("example.fastq", "fastq"))
>>> len(record)
25
>>> print(record.seq)
CCCTTCTTGTCTTCAGCGTTTCTCC
>>> print(record.letter_annotations["phred_quality"])
[26, 26, 18, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 22, 26, 26, 26, 26, 26, 26, 26, 23, 23]
假设这是 Roche 454 数据,并且从其他信息来看,您认为 TTT
应该只是 TT
。我们可以通过首先在“额外”的第三个 T
之前和之后对 SeqRecord
进行切片来创建一个新的编辑后的记录
>>> left = record[:20]
>>> print(left.seq)
CCCTTCTTGTCTTCAGCGTT
>>> print(left.letter_annotations["phred_quality"])
[26, 26, 18, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 22, 26, 26, 26, 26]
>>> right = record[21:]
>>> print(right.seq)
CTCC
>>> print(right.letter_annotations["phred_quality"])
[26, 26, 23, 23]
现在将这两个部分加在一起
>>> edited = left + right
>>> len(edited)
24
>>> print(edited.seq)
CCCTTCTTGTCTTCAGCGTTCTCC
>>> print(edited.letter_annotations["phred_quality"])
[26, 26, 18, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 22, 26, 26, 26, 26, 26, 26, 23, 23]
简单直观?我们希望如此!您可以将其简化为
>>> edited = record[:20] + record[21:]
现在,对于带有特征的示例,我们将使用 GenBank 文件。假设您有一个环状基因组
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])
>>> len(record)
9609
>>> len(record.features)
41
>>> record.dbxrefs
['Project:58037']
>>> record.annotations.keys()
dict_keys(['molecule_type', 'topology', 'data_file_division', 'date', 'accessions', 'sequence_version', 'gi', 'keywords', 'source', 'organism', 'taxonomy', 'references', 'comment'])
您可以像这样移动起点
>>> shifted = record[2000:] + record[:2000]
>>> shifted
SeqRecord(seq=Seq('GATACGCAGTCATATTTTTTACACAATTCTCTAATCCCGACAAGGTCGTAGGTC...GGA'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])
>>> len(shifted)
9609
请注意,这并不完美,因为一些注释(如数据库交叉引用、所有注释(除了分子类型)和一个特征(源特征))已经丢失
>>> len(shifted.features)
40
>>> shifted.dbxrefs
[]
>>> shifted.annotations.keys()
dict_keys(['molecule_type'])
这是因为 SeqRecord
切片步骤对保留哪些注释非常谨慎(错误地传播注释会导致重大问题)。如果您想保留数据库交叉引用或 annotations 字典,则必须显式执行此操作
>>> shifted.dbxrefs = record.dbxrefs[:]
>>> shifted.annotations = record.annotations.copy()
>>> shifted.dbxrefs
['Project:58037']
>>> shifted.annotations.keys()
dict_keys(['molecule_type', 'topology', 'data_file_division', 'date', 'accessions', 'sequence_version', 'gi', 'keywords', 'source', 'organism', 'taxonomy', 'references', 'comment'])
另请注意,在这样的示例中,您可能应该更改记录标识符,因为 NCBI 引用指的是原始未修改的序列。
反向互补 SeqRecord 对象
Biopython 1.57 中的新功能之一是 SeqRecord
对象的 reverse_complement
方法。它试图在易用性和对反向互补记录中的注释如何处理的担忧之间取得平衡。
对于序列,它使用 Seq 对象的反向互补方法。任何特征都将被转移,位置和链将被重新计算。同样,任何每字母注释也会被复制但被反转(这对典型的示例(如质量分数)是有意义的)。然而,大多数注释的转移是有问题的。
例如,如果记录 ID 是一个登录号,那么这个登录号实际上不应该应用于反向互补序列,默认情况下转移标识符很容易导致下游分析中出现微妙的数据损坏。因此,默认情况下,SeqRecord
的 id、名称、描述、注释和数据库交叉引用都不会被默认转移。
SeqRecord
对象的 reverse_complement
方法接受许多与记录属性相对应的可选参数。将这些参数设置为 True
表示复制旧值,而 False
表示删除旧值并使用默认值。您也可以提供新的所需值。
考虑这个示例记录
>>> from Bio import SeqIO
>>> rec = SeqIO.read("NC_005816.gb", "genbank")
>>> print(rec.id, len(rec), len(rec.features), len(rec.dbxrefs), len(rec.annotations))
NC_005816.1 9609 41 1 13
这里我们进行反向互补并指定一个新的标识符 - 但请注意,大多数注释都被删除了(但特征除外)
>>> rc = rec.reverse_complement(id="TESTING")
>>> print(rc.id, len(rc), len(rc.features), len(rc.dbxrefs), len(rc.annotations))
TESTING 9609 41 0 0