顺序注释对象
章 序列对象 介绍了序列类。在“上方”, Seq
类是序列记录或 SeqRecord
类,在 Bio.SeqRecord
module.此类允许更高级别的功能,例如标识符和功能(如 SeqFeature
对象)与序列关联,并在整个序列输入/输出界面中使用 Bio.SeqIO
在第四章中详细描述 序列输入/输出 .
如果您只打算使用FASTA文件等简单数据,那么您可能可以暂时跳过本章。另一方面,如果您要使用经过丰富注释的序列数据,例如来自基因库或MBE文件的序列数据,那么该信息就非常重要。
虽然本章应该涵盖与 SeqRecord
和 SeqFeature
本章中的对象,您可能还想阅读 SeqRecord
维基页面(http://biopython.org/wiki/SeqRecord)和内置文档 (Bio.SeqRecord
和 Bio.SeqFeature
):
>>> from Bio.SeqRecord import SeqRecord
>>> help(SeqRecord)
SeqRecord对象
的 SeqRecord
(序列记录)类在 Bio.SeqRecord
module.此类允许将标识符和要素等更高级的要素与序列关联(请参阅第章 序列对象 ),并且是 Bio.SeqIO
序列输入/输出接口(请参阅第章 序列输入/输出 ).
的 SeqRecord
类本身非常简单,并提供以下信息作为属性:
- .seq
序列本身,通常是
Seq
object.- .id
用于标识序列的主ID-字符串。在大多数情况下,这是一个类似于登录号的东西。
- .姓名
序列的“公共”名称/id-字符串。在某些情况下,这将与登录号相同,但它也可以是克隆名称。我认为这类似于基因库记录中的LOCUS id。
- 。描述
序列的人类可读描述或表达性名称-字符串。
- .letter_annotations
使用有关序列中字母的其他信息的(限制)字典保存按字母的注释。键是信息的名称,信息作为Python序列(即列表、数组或字符串)包含在值中,长度与序列本身相同。这通常用于质量评分(例如部分 FASTQ文件的简单质量过滤 )或二级结构信息(例如来自Stockholm/PFAM对齐文件)。
- 。注释
有关序列的附加信息的字典。键是信息的名称,信息包含在值中。这允许向序列添加更多“非结构化”信息。
- 。特点
的列表
SeqFeature
具有有关序列特征的更多结构化信息(例如基因组上基因的位置或蛋白质序列上的结构域)的对象。序列特征的结构在下面的部分中描述 特征、位置和位置对象 .- .dbxrefs
作为字符串形式的数据库交叉引用列表。
创建SeqRecord
使用 SeqRecord
对象并不是很复杂,因为所有信息都作为类的属性呈现。通常您不会创建 SeqRecord
“手工”,但改为使用 Bio.SeqIO
为您读取序列文件(请参阅第章 序列输入/输出 以及下面的例子)。然而,创造 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文件,其中包含整个序列 Yersinia pestis biovar Microtus str.91001 pCHP 1,最初从NCBI下载。该文件包含在GeneBank文件夹下或在线的Biopython单元测试中 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文件有所不同,因此通常这是不可能的。
基因库文件中的SeqRecord对象
与前面的例子一样,我们将查看整个序列 Yersinia pestis biovar Microtus str.91001 pCHP 1,最初从NCBI下载,但这次是作为基因库文件。同样,该文件包含在Gene文件夹下或在线的Biopython单元测试中 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')
的 name
来自LOCUS系列,而 id
包括版本后缀。描述来自定义行:
>>> record.id
'NC_005816.1'
>>> record.name
'NC_005816'
>>> record.description
'Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
SEN文件没有任何按字母排列的注释:
>>> record.letter_annotations
{}
大多数注释信息都记录在 annotations
词典,例如:
>>> len(record.annotations)
13
>>> record.annotations["source"]
'Yersinia pestis biovar Microtus str. 91001'
的 dbxrefs
列表将从任何PROCEMENT或DBLINK行填充:
>>> record.dbxrefs
['Project:58037']
最后,也许最有趣的是,特征表中的所有条目(例如基因或CDS特征)都被记录为 SeqFeature
中的对象 features
名单
>>> len(record.features)
41
我们再谈 SeqFeature
接下来,在部分中 特征、位置和位置对象 .
特征、位置和位置对象
SeqPerformance对象
序列特征是描述序列的重要组成部分。一旦您超越了序列本身,您需要某种方法来组织并轻松获取有关序列的更“抽象”信息。虽然可能不可能开发涵盖所有内容的通用序列要素类,但Biopython SeqFeature
类尝试封装尽可能多的有关序列的信息。该设计在很大程度上基于ESB/MBE功能表,因此如果您了解它们的外观,您可能会更容易掌握Biopython类的结构。
关于每个的关键想法 SeqFeature
目标是描述父序列上的区域,通常是 SeqRecord
object.该区域用位置对象来描述,通常是两个位置之间的范围(请参阅第节 位置和场所 下面)。
的 SeqFeature
类有许多属性,因此首先我们将列出它们及其一般功能,然后在本章的后面通过示例来展示这如何应用于现实生活中的示例。SeqPerformance的属性包括:
- .类型
这是特征类型的文本描述(例如,这类似于“CDS”或“基因”)。
- 。地点
的位置
SeqFeature
关于您正在处理的序列,请参阅部分 位置和场所 下面的SeqFeature
将其大部分功能委托给位置对象,并包括位置属性的许多快捷属性:- .ref
的简写
.location.ref
- 该位置所引用的任何(不同)参考序列。通常只是没有。- .ref_db
的简写
.location.ref_db
- 指定数据库中的任何标识符.ref
指的是。通常只是没有。- .strand
的简写
.location.strand
- 对于双链核苷酸序列,这可以是 \(1\) 对于顶部的链, \(-1\) 对于底部链, \(0\) 如果该链很重要但未知,或者None
如果不重要的话。对于蛋白质或单链序列,这是无。
- 。资格赛
这是有关该功能的其他信息的Python词典。关键是对值中包含的信息的某种简短的一字描述,而值就是实际信息。例如,限定符的公共密钥可能是“证据”,而值可能是“计算性(非实验性)”。这只是一种让查看该功能的人知道它不是实验性的(i. e.在潮湿的实验室中)确认。请注意,其他值将是字符串列表(即使只有一个字符串)。这是ESB/MBE文件中特征表的反映。
- .sub_features
这曾经用于表示具有复杂位置的特征,例如ESB/MBE文件中的“连接”。随着
CompoundLocation
对象,现在应该忽略。
位置和场所
关于每个的关键想法 SeqFeature
对象是描述父序列上的区域,为此我们使用位置对象,通常描述两个位置之间的范围。两个人试图澄清我们正在使用的术语:
- 位置
这指的是序列上的单个位置,该位置可能是模糊的,也可能不是模糊的。例如,5、20、
<100
和>200
都是位置。- 位置
位置是以某些位置为界的序列区域。例如
5..20
(i. e.(5 - 20)是一个位置。
我之所以提到这一点,是因为有时我会混淆这两者。
SimpleLocate对象
除非你使用真核基因,否则大多数 SeqFeature
位置非常简单-您只需要开始和结束坐标和一条线。这本质上是所有基本的 SimpleLocation
对象确实如此。
当然,在实践中,事情可能会更加复杂。首先,我们必须处理由多个地区组成的大院位置。其次,位置本身可能是模糊的(不准确的)。
CompoundLocate对象
Biopython 1.62引入了 CompoundLocation
作为重组由多个地区组成的复杂地点的代表方式的一部分。主要用途是处理EmbL/Gene文件中的“加入”位置。
模糊位置
到目前为止,我们只使用了简单的位置。处理特征位置的一个复杂之处在于位置本身。在生物学中,很多时候事情并不完全确定(就像我们湿实验室生物学家试图使它们确定一样!)。例如,你可能做了一个二核苷酸引发实验,发现mRNA转录的起始点是两个位点之一。这是非常有用的信息,但复杂的是如何将其表示为一个位置。为了帮助我们解决这个问题,我们有了模糊位置的概念。基本上有几种类型的模糊位置,所以我们有五个类别来处理它们:
- ExactPosition
顾名思义,这个类代表沿着序列指定为精确的位置。这只是一个数字,您可以通过查看
position
对象的属性。- BeforePosition
这个类表示出现在某个指定位置之前的模糊位置。在GenBank/EMBL标记法中,这表示为类似于
<13
表示实际位置位于小于13的某个位置。若要获取指定的上边界,请查看position
对象的属性。- AfterPosition
违反
BeforePosition
,此类表示出现在某个指定站点之后的位置。这在基因库中表示为>13
,和喜欢BeforePosition
,您可以通过查看position
对象的属性。- WithinPosition
此类有时用于基因库/MBE位置,建模出现在两个指定核苷酸之间某处的位置。在SEN/MBE符号中,这将被表示为
(1.5)
,以表示该位置位于1到5范围内的某个地方。- OneOfPosition
偶尔用于基因库/MBE位置,此类处理存在多个可能值的位置,例如,如果起始密码子不清楚并且那里有两个候选基因起始,则可以使用它。或者,这可能会被明确地视为两个相关的基因特征。
- UnknownPosition
此类处理未知位置的位置。这不在ESB/MBE中使用,但对应于“?”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)]
我们可以使用位置的开始和结束属性访问模糊开始和结束位置:
>>> 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']
请注意,用连接定义的来自基因库或MBE文件的基因和CDS特征是exon的联合-它们不覆盖任何intron。
由特征或位置描述的序列
A SeqFeature
或位置对象不直接包含序列,而是包含位置(请参阅第节 位置和场所 )描述了如何从父序列中获取该内容。例如,考虑具有位置的(短)基因序列 5:18
在反向链上,在使用基于1的计数的基因库/MBE标记法中将是 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
长度是组成区域的总和。
比较
的 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的旧版本会使用Python的默认对象比较 SeqRecord
,意思是 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
对象,以便它可以指定引用引用所引用的序列上的特定位置。例如,您可能有一本期刊,该期刊正在处理位于BAT上的特定基因,并希望指定它仅准确引用该位置。的 location
是一个潜在的模糊位置,如部分所述 位置和场所 .
任何引用对象都作为列表存储在 SeqRecord
对象的 annotations
字典在关键“参考”下。这就是全部。引用旨在易于处理,并且希望足够通用以涵盖大量用例。
format方法
的 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
这 format
方法采用单个强制参数,即由支持的大小写字符串 Bio.SeqIO
作为输出格式(请参阅第章 序列输入/输出 ).然而,某些文件格式 Bio.SeqIO
可以写信给 require 不止一条记录(通常是多个序列比对格式的情况),因此无法通过此工作 format()
法另见章节 将SeqRecord对象作为格式化字符串 .
切片SeqRecord
你可以切片 SeqRecord
,给你一个新的 SeqRecord
仅涵盖序列的一部分。这里重要的是,任何按字母计算的注释也会被切片,并且任何完全属于新序列的特征都会被保留(调整其位置)。
例如,采用之前使用的相同的SEN文件:
>>> 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
.如果您直接查看基因库文件,您会发现该基因/CDS具有位置字符串 4343..4780
,或在Python中计数 4342:4780
.通过查看文件,您可以发现这些是文件中的第十二和第十三个条目,因此在Python零基计数中,它们是条目 \(11\) 和 \(12\) 在 features
列表:
>>> 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
我们的子记录只有两个特征,基因和CDS条目 YP_pPCP05
:
>>> 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
[]
您可能希望保留其他条目,例如生物体?小心复制整个注释字典,因为在这种情况下,您的部分序列不再是环形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'
这很好地说明了问题,我们的新子记录是 not 该载体的完整序列,所以描述是错误的!让我们修复这个问题,然后使用 format
上述第节中描述的方法 format方法 :
>>> 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...
添加SeqRecord对象
您可以添加 SeqRecord
对象放在一起,给出新的 SeqRecord
.这里重要的是,还添加了任何常见的按字母注释,所有功能都被保留(调整了其位置),并且还保留了任何其他常见注释(例如id、名称和描述)。
对于具有按字母注释的示例,我们将使用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]
假设这是罗氏454数据,并且从其他信息中您认为 TTT
应该只有 TT
.我们可以创建一个新的编辑过的记录, SeqRecord
“额外”三分之一之前和之后 T
:
>>> 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:]
现在,对于具有功能的示例,我们将使用SEN文件。假设你有一个环形基因组:
>>> 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
切片步骤对于保留的注释很谨慎(错误地传播注释可能会导致重大问题)。如果您想保留数据库交叉引用或注释字典,则必须显式地完成:
>>> 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引用引用指的是 original 未修改的序列。
反向补充SeqRecord对象
Biopython 1.57中的新功能之一是 SeqRecord
对象的 reverse_complement
法这试图平衡易用性与对反向补记录中的注释处理的担忧。
对于序列,这使用Seq对象的反向补数方法。任何要素都将随着重新计算的位置和股而转移。同样,任何按字母计算的注释也会被复制但颠倒(这对于质量分数等典型示例是有意义的)。然而,大多数注释的传输都存在问题。
例如,如果记录ID是一个登录项,则该登录项不应真正适用于反向补序列,并且默认传输标识符可能很容易导致下游分析中微妙的数据损坏。因此,默认情况下, SeqRecord
的id、名称、描述、注释和数据库交叉引用全部 not 默认转移。
的 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