Biopython 序列注释对象实战：SeqRecord、SeqFeature 与 Location 完整指南

科研【免费下载链接】biopythonOfficial git repository for Biopython (originally converted from CVS)项目地址https://gitcode.com/gh_mirrors/bi/biopython点击查看免费下载导读本文围绕 Biopython 的序列注释体系展开核心是「序列之上的一层」——Bio.SeqRecord.SeqRecord与Bio.SeqFeature.SeqFeature两类对象。它们允许把标识符id/name/description、特征表features、逐字母注释letter_annotations、杂项注释annotations和数据库交叉引用dbxrefs与序列本身关联起来是Bio.SeqIO读写 GenBank、EMBL、FASTA、FASTQ 等格式时的基本数据类型。读完本文你将掌握如何手工或从文件构造SeqRecord如何理解SeqFeature、SimpleLocation、CompoundLocation与六类模糊位置fuzzy position对象以及format()、切片、相加、反向互补等操作对注释的保留与调整规则——这些知识是处理带注释序列数据如 GenBank 基因组记录的必备技能。本文对应的原始章节为 Doc/Tutorial/chapter_seq_annot.rst核心实现位于 Bio/SeqRecord.py 与 Bio/SeqFeature.py测试用例可参考 Tests/test_SeqRecord.py。SeqRecord 对象序列与注释的容器SeqRecord类定义在Bio.SeqRecord模块中位于Seq类「之上」一层。它允许将标识符和特征等高层信息与序列关联并贯穿于Bio.SeqIO序列读写接口详见 Doc/Tutorial/chapter_seqio.rst。如果只处理 FASTA 这类简单数据可以暂时跳过本章但若使用 GenBank 或 EMBL 等富含注释的数据这一章的内容就相当重要。SeqRecord本身相当简单信息全部以属性形式暴露属性类型说明.seqSeq对象序列本身.id字符串主标识符通常是 accession 号.name字符串「通用」名称/id可能与 accession 相同也可能是克隆名可类比 GenBank 记录的 LOCUS id.description字符串可读的描述或表达性名称.letter_annotations受限字典逐字母注释键为信息名称值为与序列等长的 Python 序列list/tuple/str常用于质量分数、二级结构等.annotations字典关于整条序列的杂项注释.featuresSeqFeature列表结构化的特征信息如基因位置、蛋白质结构域.dbxrefs字符串列表数据库交叉引用其中letter_annotations使用私有类_RestrictedDict见 Bio/SeqRecord.py它强制值必须是 Python 序列且长度必须与序列长度一致否则抛出TypeError。例如对长度 5 的序列写入hello world会报错Any per-letter annotation should be a Python sequence (list, tuple or string) of the same length as the biological sequence, here 5.这一设计可在早期拦截「注释与序列长度不匹配」这类常见错误。创建 SeqRecord从零开始创建创建SeqRecord至少要有一个Seq对象 from Bio.Seq import Seq simple_seq Seq(GATC) from Bio.SeqRecord import SeqRecord simple_seq_r SeqRecord(simple_seq)也可以把 id、name、description 直接传给初始化函数不传时会被设为表示「未知」的字符串之后可再修改 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 simple_seq_r.seq Seq(GATC)从源码看SeqRecord.__init__的完整签名Bio/SeqRecord.py支持seq必填Seq或MutableSeq、id、name、description、dbxrefs、features、annotations、letter_annotations全部参数。如果要输出到文件标识符非常重要建议创建时直接带上 simple_seq_r SeqRecord(simple_seq, idAC12345)添加杂项注释直接操作annotations字典即可 simple_seq_r.annotations[evidence] None. I just made it up. print(simple_seq_r.annotations) {evidence: None. I just made it up.}逐字母注释的操作类似值可以是与序列等长的字符串、列表或元组 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 文件创建以Yersinia pestis biovar Microtusstr. 91001 质粒 pPCP1 的全序列 FASTA 文件为例仓库内置在 Tests/GenBank/NC_005816.fna。文件开头是gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus ... pPCP1, complete sequence TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGGGGGTAATCTGCTCTCC ...对于只含单条记录的文件用Bio.SeqIO.read(...)而不是循环遍历的parse(...) from Bio import SeqIO record SeqIO.read(NC_005816.fna, fasta) record SeqRecord(seqSeq(TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG), idgi|45478711|ref|NC_005816.1|, namegi|45478711|ref|NC_005816.1|, descriptiongi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence, dbxrefs[])逐个查看关键属性 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整行标题用作description。这是有意为之部分出于向后兼容并且对有多个词标题的 FASTA 文件也说得通。其余注释属性在读取 FASTA 时都不会被填充 record.dbxrefs [] record.annotations {} record.letter_annotations {} record.features []文档也提醒虽然 NCBI 对 FASTA 标题行有一套约定可从中解析 GI 号、accession 等但其他来源的 FASTA 格式各异因此这种解析在一般情况下不可行。从 GenBank 文件创建同样使用 NC_005816 质粒但这次是 GenBank 格式仓库内置在 Tests/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 ...读取代码与 FASTA 几乎一致仅格式名不同 record SeqIO.read(NC_005816.gb, genbank) record SeqRecord(seqSeq(TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG), idNC_005816.1, nameNC_005816, descriptionYersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence, dbxrefs[Project:58037])GenBank 与 FASTA 的解析约定不同name来自 LOCUS 行id带版本后缀description来自 DEFINITION 行 record.id NC_005816.1 record.name NC_005816 record.description Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequenceGenBank 文件没有逐字母注释大部分注释信息进入annotations字典 record.letter_annotations {} len(record.annotations) 13 record.annotations[source] Yersinia pestis biovar Microtus str. 91001dbxrefs列表从 PROJECT / DBLINK 行填充 record.dbxrefs [Project:58037]最有趣的是特征表所有 entries基因、CDS 等都成为features列表中的SeqFeature对象 len(record.features) 41Feature、Location 与 Position 对象SeqFeature 对象序列特征是描述序列不可或缺的部分。SeqFeature的设计大量借鉴了 GenBank/EMBL 特征表如果你熟悉它们的结构理解 Biopython 类会更容易。每个SeqFeature的核心思想是描述父序列通常是SeqRecord上的一个区域该区域由一个 location 对象通常是两个位置之间的范围描述。SeqFeature的属性包括.type特征的文本类型描述如CDS、gene。.location特征在序列上的位置见下文「Positions and locations」。SeqFeature将大量功能委托给 location 对象并提供几个快捷属性.reflocation.ref的简写——location 所引用的不同的参考序列通常为None。.ref_dblocation.ref_db的简写——ref中标识符所属的数据库通常为None。.strandlocation.strand的简写——特征所在链。双链核酸序列上1表示顶链-1表示底链0表示链重要但未知None表示无关蛋白质或单链序列即为None。.qualifiers特征的附加信息字典。键是简短的词值是具体信息。注意值总是字符串列表即使只有一个字符串这反映了 GenBank/EMBL 特征表的结构。.sub_features曾用于表示joins等复杂位置随着CompoundLocation对象的引入已废弃现在应忽略。从源码看Bio/SeqFeature.pySeqFeature的__init__接受location必须是SimpleLocation、CompoundLocation或None、type、id、qualifiers传入sub_features会直接抛出TypeError提示改用CompoundLocation。其__str__会输出type:、location:和按字母排序的qualifiers键值对这正是打印特征时看到的格式。Positions 和 locations 术语position位置序列上的单个位置可能是模糊fuzzy的也可能不是。例如 5、20、100、200都是位置。location位置区域由若干位置界定的序列区域。例如5..20就是一个 location。SimpleLocation 对象除非处理真核基因大多数SeqFeature的 location 都极其简单——只需要起始坐标、结束坐标和链。SimpleLocation基本就做这件事。注意编号遵循 Python 的规则GenBank 中 1-based 的123..150对应零基的[122:150]参见 Bio/SeqFeature.py 的 docstring。SimpleLocation的start与end支持传入整数自动转为ExactPosition或任意*Position对象也可附带strand、ref、ref_dbBio/SeqFeature.py。对于长度为 n 的父序列start与end必须满足0 start end n即使是反链特征也把start理解为「最左最小边界」end理解为「最右最大边界」这在复合位置中尤其重要。CompoundLocation 对象Biopython 1.62 引入CompoundLocation用于表示由多个区域组成的复杂位置主要场景是 EMBL/GenBank 文件中的join位置。它替代了旧的sub_features机制。源码实现Bio/SeqFeature.py要求至少两个SimpleLocation部件并支持join、order等 operator整体strand由各部件自动推导——若所有部件链一致则返回该值混合链时返回None。CompoundLocation还支持用运算符把多个SimpleLocation直观地拼接起来构造。模糊位置Fuzzy Positions生物学中很多位置并非完全确定例如双核苷酸引物实验发现 mRNA 转录起始位点有两个可能位点。Biopython 用以下位置类处理不确定性类含义GenBank/EMBL 表示ExactPosition精确定位读position属性直接数字BeforePosition真实位置在指定点之前13AfterPosition真实位置在指定点之后13WithinPosition位置在某两个核苷酸之间(1.5)OneOfPosition位置有多个可能取值如起始密码子不清—UnknownPosition位置未知对应 UniProt 的?非 GenBank/EMBL 使用源码中还有UncertainPosition类见 Bio/SeqFeature.py。创建带模糊端点的 location 示例 from Bio import SeqFeature start_pos SeqFeature.AfterPosition(5) end_pos SeqFeature.BetweenPosition(9, left8, right9) my_location SeqFeature.SimpleLocation(start_pos, end_pos)注意自 Biopython 1.59 起BetweenPosition和WithinPosition必须显式指定切片等操作应使用的整数位置——起始位置一般取较小left值结束位置一般取较大right值。打印 location 可以得到不错的表示 print(my_location) [5:(8^9)]通过start/end属性访问模糊端点 my_location.start AfterPosition(5) print(my_location.start) 5 my_location.end BetweenPosition(9, left8, right9) print(my_location.end) (8^9)位置类是整数的子类可以直接转成整数用 int(my_location.start) 5 int(my_location.end) 9如果想避开模糊位置的细节直接给整数即可得到ExactPosition exact_location SeqFeature.SimpleLocation(5, 9) print(exact_location) [5:9] exact_location.start ExactPosition(5) int(exact_location.start) 5Location 测试判断某位置是否落在特征内可以用 Python 的in关键字对SeqFeature或 location 对象做成员测试判断父序列坐标上的某个碱基/残基是否在该特征内。例如有个 SNP 位于索引 4350Python 计数暴力遍历所有特征 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 文件中用 joins 定义的 gene 与 CDS 特征是各外显子的并集不含内含子。源码中SeqFeature.__contains__Bio/SeqFeature.py正是委托给self.location的成员测试对CompoundLocation而言部件之间的空隙不会被计入。模糊位置如BeforePosition也可能需要额外小心测试用例展示了相应行为。特征或 location 描述的序列extract()SeqFeature或 location 对象本身不直接包含序列location 描述的是如何从父序列中取到它。例如一个反链基因location 为5:18在 GenBank/EMBL 的 1-based 记法中是complement(6..18) from Bio.Seq import Seq from Bio.SeqFeature import SeqFeature, SimpleLocation seq Seq(ACCGAGACGGCAAAGGCTAGCATAGGTATGAGACTTCCTTCCTGCCAGTGCTGAGGAACTGGGAGCCTAC) feature SeqFeature(SimpleLocation(5, 18, strand-1), typegene)可以手动切片再取反向互补——location 的 start/end 是整数似的可以直接用 feature_seq seq[feature.location.start : feature.location.end].reverse_complement() print(feature_seq) AGCCTTTGCCGTC但对于复合特征joins手动处理很麻烦。SeqFeature提供extract方法Bio/SeqFeature.py一步到位自 Biopython 1.78 起还可通过传入参考序列字典支持反式剪接trans-splicing feature_seq feature.extract(seq) print(feature_seq) AGCCTTTGCCGTCSeqFeature或 location 的长度与其描述的序列区域一致 print(len(feature_seq)) 13 print(len(feature)) 13 print(len(feature.location)) 13对SimpleLocation长度就是 start 与 end 之差对CompoundLocation长度是各组成部分之和内含子等空隙不计入保证len(f)与len(f.extract(parent_seq))一致。若 location 为None解析到非法位置时可能出现extract()会抛出ValueError提示检查序列文件的 location。SeqRecord 比较刻意不实现两个看似相同的记录能否直接比较例如 from Bio.Seq import Seq from Bio.SeqRecord import SeqRecord record1 SeqRecord(Seq(ACGT), idtest) record2 SeqRecord(Seq(ACGT), idtest)旧版 Biopython 使用 Python 默认的对象比较record1 record2只在指向同一内存对象时为True上例会返回False很容易让人踩坑。自 Biopython 1.67 起SeqRecord的比较会显式抛出错误 record1 record2 Traceback (most recent call last): ... NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest.从源码看Bio/SeqRecord.py 中__eq__、__ne__、__lt__、__le__、__gt__、__ge__全部抛出NotImplementedError消息字符串定义在 Bio/SeqRecord.py。正确的做法是逐项比较你关心的属性 record1.id record2.id True record1.seq record2.seq True注意复杂对象的比较很快就会变得麻烦可参考 Doc/Tutorial/chapter_seq_objects.rst 中的序列比较一节。Tests/test_SeqRecord.py中的test_eq_exception、test_lt_exception等测试用例即覆盖这些行为。References参考文献对象Bio.SeqFeature.Reference类以对象属性的方式存储期刊文献等引用信息属性包括journal、title、authors、medline_id、pubmed_id、comment以及一个location对象用于指定该文献所对应的序列区域可以是模糊位置。源码见 Bio/SeqFeature.py其中location是Location对象列表为空则默认覆盖整条序列。所有引用对象以列表形式存放在SeqRecord.annotations字典的references键下。References 的设计初衷就是「容易处理、足够通用」。format() 方法把 SeqRecord 变成文件格式字符串SeqRecord.format()返回将记录按Bio.SeqIO支持的输出格式如 FASTA格式化后的字符串 from Bio.Seq import Seq from Bio.SeqRecord import SeqRecord record SeqRecord( ... Seq( ... MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD ... GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK ... NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM ... SSAC ... ), ... idgi|14150838|gb|AAK54648.1|AF376133_1, ... descriptionchalcone synthase [Cucumis sativus], ... ) print(record.format(fasta))输出gi|14150838|gb|AAK54648.1|AF376133_1 chalcone synthase [Cucumis sativus] MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM SSACformat()接受一个必填参数Bio.SeqIO支持的小写输出格式名详见 Doc/Tutorial/chapter_seqio.rst。但有些Bio.SeqIO能写出的格式需要多条记录典型如多序列比对格式因此无法通过单个SeqRecord的format()使用。从源码看Bio/SeqRecord.pyformat()委托给__format__后者通过SeqIO._FormatToWriter[format_spec].to_string(self)实现二进制格式会抛出ValueError。__format__还支持format(record, fasta)与 f-string如f{record:fasta}两种写法。切片 SeqRecord对SeqRecord切片会得到覆盖部分序列的新SeqRecord。关键规则是逐字母注释同步切片完全落在新序列内的特征保留且位置已调整。继续使用之前的 GenBank 文件 from Bio import SeqIO record SeqIO.read(NC_005816.gb, genbank) len(record) 9609 len(record.features) 41聚焦pim基因YP_pPCP05。在 GenBank 文件中该 gene/CDS 的 location 为4343..4780即 Python 计数4342:4780从文件可数出它们是第 12、13 个条目即features列表的索引 11 与 12这里文档演示用索引 20 与 21 打印实际索引取决于文件具体内容 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: product, Value: [pesticin immunity protein] Key: protein_id, Value: [NP_995571.1] Key: transl_table, Value: [11] ...现在切片父记录 4300 到 4800覆盖pimgene/CDS看看保留了多少特征 sub_record record[4300:4800] len(sub_record) 500 len(sub_record.features) 2子记录只剩两个特征——YP_pPCP05的 gene 与 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 ...源码实现Bio/SeqRecord.py确认了这些行为name、id、description原样保留逐字母注释随子序列切片不使用步长时完全落入子序列的特征被保留并调整位置_shift除 molecule_type 外annotations字典与dbxrefs列表不会带过去因为对子序列来说它们可能已不再适用。molecule_type之所以例外是因为输出 GenBank/EMBL 时仍然需要它 sub_record.annotations {molecule_type: DNA} sub_record.dbxrefs []是否要保留 organism 等其他条目由你自己决定。切忌整体复制 annotations 字典比如本例中部分序列已不再是环状 DNA而是线性的应手动修正 sub_record.annotations[topology] linearid、name、description出于实用考虑被保留但这也带来问题新子记录并不是质粒的完整序列原 description 已不正确。修改后再用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 场景下切片会同步裁剪逐字母注释质量分数相关示例可参考 Doc/Tutorial/chapter_seqio.rst 中的 FASTQ 剪除引物/接头章节。测试用例见Tests/test_SeqRecord.py的test_slicing、test_slice_variants、test_slice_simple、test_slice_zero等。相加 SeqRecordSeqRecord支持相加得到新的SeqRecord。规则是共同的逐字母注释也会相加所有特征保留位置已调整其他共同注释id、name、description保留。先看逐字母注释的例子用 FASTQ 文件的第一条记录测试数据见 Tests/Quality/example.fastq 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 前后各切片一部分 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:]。再看带特征的例子。假设有一个环状基因组Tests/GenBank/NC_005816.gb可以这样移动原点 from Bio import SeqIO record SeqIO.read(NC_005816.gb, genbank) 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] len(shifted) 9609但这并不完美数据库交叉引用、除 molecule_type 外的全部注释、以及 source 特征都丢失了 len(shifted.features) 40 shifted.dbxrefs [] shifted.annotations.keys() dict_keys([molecule_type])原因是SeqRecord切片在保留注释时非常谨慎错误地传播注释可能造成大问题。若想保留 dbxrefs 或 annotations必须显式拷贝 shifted.dbxrefs record.dbxrefs[:] shifted.annotations record.annotations.copy() shifted.dbxrefs [Project:58037]另外这类操作后最好修改记录标识符因为 NCBI 引用指向的是原始的未修改序列。源码实现Bio/SeqRecord.py显示相加两个SeqRecord时第二个记录的每个特征通过f._shift(length)平移到新坐标只有两者相同的 id/name/description 才保留annotations只保留两边一致的键逐字母注释若两边都有同名键则拼接长度校验会被绕过以避免中间态误报。若与普通字符串/Seq相加则逐字母注释会丢失见 Bio/SeqRecord.py。测试用例见Tests/test_SeqRecord.py的test_add_seqrecord、test_add_seq、test_slice_add_shift。反向互补 SeqRecordreverse_complement()是 Biopython 1.57 加入的功能在易用性与「注释怎么办」的担忧之间做了平衡序列使用Seq对象的反向互补方法所有特征转移location 与 strand 重算逐字母注释复制并反转对质量分数这类例子是合理的但绝大多数注释的转移有问题——例如若记录 id 是 accession它并不真正适用于反向互补后的序列默认转移标识符可能在下游分析中造成隐蔽的数据损坏。因此默认不转移id、name、description、annotations和 dbxrefs。该方法接受对应属性的可选参数设为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(idTESTING) print(rc.id, len(rc), len(rc.features), len(rc.dbxrefs), len(rc.annotations)) TESTING 9609 41 0 0从源码看Bio/SeqRecord.py默认签名为reverse_complement(idFalse, nameFalse, descriptionFalse, featuresTrue, annotationsFalse, letter_annotationsTrue, dbxrefsFalse)——只有 features 与 letter_annotations 默认保留。逐字母注释反转后Tests/Quality/solexa_faked.fastq的示例中 Solexa 质量值会由[40, 39, ..., -5]变为[-5, -4, ..., 40]。特征在反向互补后 strand 翻转、位置镜像例如 pBAD30 质粒Tests/GenBank/pBAD30.gb的第一个 CDS 从1081:1960变为2963:3842由len(plasmid) - 1081 3842、len(plasmid) - 1960 2963可验证。需要留意两点若特征注释包含链特异信息如 SNP 的碱基变化这些信息不会被修正反向互补后需自行处理对蛋白质记录做反向互补会抛出异常例如molecule_type为protein时。测试用例见Tests/test_SeqRecord.py的test_reverse_complement_seq、test_reverse_complement_mutable_seq。总结与实战要点以本仓库为例把整套流程串起来就是用Bio.SeqIO读取 Tests/GenBank/NC_005816.gb 得到SeqRecord→ 遍历record.features用in测试 SNP 位置 → 用feature.extract(seq)提取基因序列 → 用切片与相加移动环状基因组原点 → 用reverse_complement()生成反链记录 → 用format(genbank)/format(fasta)输出。核心要点可归纳为注释分五层id/name/description字符串、annotations杂项字典、letter_annotations受限等长字典、featuresSeqFeature列表、dbxrefs字符串列表。位置体系分层SeqFeature→locationSimpleLocation/CompoundLocation→position六类模糊位置类整数坐标可直接当 int 使用。操作要谨慎切片与相加默认丢弃annotations/dbxrefs仅保留 molecule_type需要时显式copy()SeqRecord比较刻意未实现反向互补默认丢弃记录级注释。优先用方法而非手写extract()处理复合与模糊位置format()借助Bio.SeqIO输出遇到二进制格式会抛出ValueError。赞分享科研【免费下载链接】biopythonOfficial git repository for Biopython (originally converted from CVS)项目地址https://gitcode.com/gh_mirrors/bi/biopython点击查看免费下载相关推荐Spring-Loaded在Grails中的应用提升开发效率的完整方案Spring Loaded在Grails中的应用提升开发效率的完整方案 Spring Loaded是一款强大的Java代理工具能够实现运行时类重新加载功能开发工具Biopython序列特征提取从基因注释到功能预测的完整指南Biopython序列特征提取从基因注释到功能预测的完整指南 Biopython是生物信息学领域最强大的Python工具包之一专门用于处理生物序列数据。在基科研Biopython 序列 Motif 分析完全指南Bio.motifs 对象模型、PWM/PSSM 与多格式解析实战Biopython 序列 Motif 分析完全指南Bio.motifs 对象模型、PWM/PSSM 与多格式解析实战 导读 Bio.motifs 是 Biop科研上一篇无TPM也能装Win11Rufus三分钟做出系统启动盘下一篇掌握知识主权用dedao-dl构建你的终身数字知识库创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考