Biopython系列教程(二)序列与注释对象详解

在上一篇《快速上手与基础入门》中，我们了解了 Biopython 的安装和基本数据结构。本篇我们将深入探讨 Biopython 中两个核心对象：Seq 和 SeqRecord。它们是处理生物序列及其注释的基础，也是后续学习序列输入输出和高级分析的前提。

生物序列是生物信息学的核心。Biopython 提供了 Seq 对象来存储和操作序列，它本质上是一个带有生物学方法的字符串。

from Bio.Seq import Seqmy_seq = Seq("GATCG")my_seq# Seq('GATCG')

Seq 对象支持大部分 Python 字符串的方法，例如获取长度、遍历、索引和切片。

# 遍历序列for index, letter in enumerate(my_seq):    print("%i %s" % (index, letter))# 0 G# ... # 获取长度print(len(my_seq))    # 5# 访问元素: 注意Python从零开始计数print(my_seq[0])    # G（第一个字母）print(my_seq[2])    # T（第三个字母）

切片操作与字符串完全一致

from Bio.Seq import Seqmy_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")my_seq[4:12]    # Seq('GATGGGCC')# 用步长提取特定位置的碱基: start, stop, stride (步长默认值为1)print(my_seq[0::3])    # GCTGTAGTAAGprint(my_seq[1::3])    # AGGCATGCATC# 使用 [::-1] 可以反转序列print(my_seq[::-1])     # CGCTAAAAGCTAGGATATATCCGGGTAGCTAG

Seq 的 .count() 方法与字符串相同，注意它是非重叠计数：

print("AAAA".count("AA"))    # 2 (非重叠)print(Seq("AAAA").count("AA"))    # 2

若要计算 GC 含量，可以手动计算：

from Bio.Seq import Seqmy_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")100 * (my_seq.count("G") + my_seq.count("C")) / len(my_seq)

但更推荐使用 Bio.SeqUtils.gc_fraction()，它能处理混合大小写序列和简并碱基：

from Bio.SeqUtils import gc_fractionmy_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")print(gc_fraction(my_seq))    # 0.46875my_seq = Seq("SATcgATgggcCTATATAGGATCGAAAATCGC")print(gc_fraction(my_seq))    # 0.46875

两个 Seq 对象可以直接用 + 拼接：

from Bio.Seq import Seqseq1 = Seq("ACGT")seq2 = Seq("AACCGG")seq1 + seq2# Seq('ACGTAACCGG')

from Bio.Seq import Seqlist_of_seqs = [Seq("ACGT"), Seq("AACC"), Seq("GGTT")]concatenated = Seq("")for s in list_of_seqs:    concatenated += sconcatenated# Seq('ACGTAACCGGTT')

from Bio.Seq import Seqcontigs = [Seq("ATG"), Seq("ATCCCG"), Seq("TTGCA")]spacer = Seq("N" * 10)spacer.join(contigs)Seq('ATGNNNNNNNNNNATCCCGNNNNNNNNNNTTGCA')

与字符串一样，Seq 有 .upper() 和 .lower() 方法，方便进行大小写不敏感匹配：

from Bio.Seq import Seqdna_seq = Seq("acgtACGT")dna_seq            # Seq('acgtACGT')dna_seq.upper()    # Seq('ACGTACGT')dna_seq.lower()    # Seq('acgtacgt')print("GTAC" in dna_seq)        # Falseprint("GTAC" in dna_seq.upper())# True

对于核苷酸序列，可以使用 .complement() 和 .reverse_complement()：

from Bio.Seq import Seqmy_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC")my_seq.complement()# Seq('CTAGCTACCCGGATATATCCTAGCTTTTAGCG')my_seq.reverse_complement()# Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC')

在生物学中，转录通常以编码链为模板，将 T 替换为 U。

实际的生物转录过程从模板链进行反向互补提供 mRNA。但在 Biopython 和一般生信中，通常直接使用编码链，因为只需通过切换 T 就可以获得 mRNA 序列，Biopython 提供了 .transcribe() 方法：

from Bio.Seq import Seqcoding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")template_dna = coding_dna.reverse_complement()coding_dna    # Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG')template_dna  # Seq('CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT')messenger_rna = coding_dna.transcribe()messenger_rna # Seq('AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG')

若你拥有模板链，需要先取反向互补再转录：

mrna = template_dna.reverse_complement().transcribe()

从 mRNA 恢复到编码链 DNA 使用 .back_transcribe()（U→T）：

from Bio.Seq import Seqmessenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG")messenger_rna.back_transcribe()# Seq('ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG')

将 mRNA 翻译成蛋白质用 .translate()，默认使用标准遗传密码：

from Bio.Seq import Seqmessenger_rna = Seq("AUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG")messenger_rna.translate()# Seq('MAIVMGR*KGAR*')

直接翻译编码链 DNA 也可以（自动转录后再翻译）：

coding_dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")coding_dna.translate()# Seq('MAIVMGR*KGAR*')

1）指定遗传密码表

Biopython 支持 NCBI 的所有遗传密码表 (https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi)。默认情况下，翻译将使用标准遗传密码 (NCBI table id 1)。

coding_dna.translate(table="Vertebrate Mitochondrial")# Seq('MAIVMGRWKGAR*')coding_dna.translate(table=2)# Seq('MAIVMGRWKGAR*')

2）翻译到第一个终止密码子

使用 to_stop=True 可以在第一个框内终止密码子处停止，并且不包含终止符号：

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

如果不喜欢默认星号，甚至可以指定停止符号：

coding_dna.translate(table=2, stop_symbol="@")# Seq('MAIVMGRWKGAR@')

Biopython 在 Bio.Data.CodonTable 中提供了所有遗传密码表，可以查看详细信息：

from Bio.Data import CodonTablestandard_table = CodonTable.unambiguous_dna_by_name["Standard"]mito_table = CodonTable.unambiguous_dna_by_name["Vertebrate Mitochondrial"]

或者，这些表分别标记有 ID 1 和 2：

from Bio.Data import CodonTablestandard_table = CodonTable.unambiguous_dna_by_id[1]mito_table = CodonTable.unambiguous_dna_by_id[2]

从 Biopython 1.65 开始，Seq 对象的比较只比较序列本身，不考虑分子类型（像 Python 字符串一样进行比较）。因此：

from Bio.Seq import Seqseq1 = Seq("ACGT")"ACGT" == seq1    # Trueseq1 == "ACGT"    # True

有时只知道序列长度而不知道具体碱基（例如 GenBank 中只提供了配置信息）。可以创建长度为 N 的未知序列：

from Bio.Seq import Sequnknown_seq = Seq(None, 10)print(len(unknown_seq)    # 10# print(unknown_seq)    # 抛出 UndefinedSequenceError

部分已知序列可以用字典指定已知片段的位置：

seq = Seq({117512683: "TTGAAAACCTGAATGTGAGAGTCAGTCAAGGATAGT"}, length=159345973)# 键是已知序列段的开始坐标，值是相应的序列内容

如果至少有一个序列部分或完全未定义，则也可以通过附加序列来创建部分定义的序列：

seq = Seq("ACGT")undefined_seq = Seq(None, length=10)seq + undefined_seq + seq# Seq({0: 'ACGT', 14: 'ACGT'}, length=18)

Seq 对象是不可变的（immutable），一旦创建就不能修改。

from Bio.Seq import Seqmy_seq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")my_seq[5] = "G"# TypeError: 'Seq' object does not support item assignment

如果需要修改序列（如模拟点突变），可以使用 MutableSeq：

from Bio.Seq import MutableSeqmutable_seq = MutableSeq(my_seq)# 或者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 不能作为字典键，因为它是可变的。

完成修改后，可以转回不可变的 Seq 对象：

new_seq = Seq(mutable_seq)

Seq 对象提供了 find、rfind、index、rindex 方法，与字符串类似，但支持多种输入类型 (str、bytes、bytearray、Seq、MutableSeq)：

from Bio.Seq import Seq, MutableSeqseq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")seq.index("ATGGGCCGC") # 9seq.find("AT")         # 3seq.index("ACTG")     # ValueError: subsection not found (未找到)seq.find("ACTG")      # -1 (未找到)seq.rindex("CC")      # 29 (从右开始查找)seq.rfind("CC")       # 29 (从右开始查找)

若想同时搜索多个子序列，可以使用 .search() 方法，它返回一个迭代器：

seq = Seq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA")for index, sub in seq.search(["CC", "GGG"]):    print(index, sub)# 1 CC# 11 GGG# 14 CC# 23 GGG# 28 CC# 29 CC

如果你更喜欢函数式编程，Bio.Seq 模块也提供了直接操作字符串的函数，它们接受字符串、Seq 或 MutableSeq 对象：

from Bio.Seq import reverse_complement, transcribe, back_transcribe, translatemy_string = "GCTGTTATGGGTCGTTGGAAGGGTGGTCGTGCTGCTGGTTAG"reverse_complement(my_string)# 'CTAACCAGCAGCACGACCACCCTTCCAACGACCCATAACAGC'transcribe(my_string)# 'GCUGUUAUGGGUCGUUGGAAGGGUGGUCGUGCUGCUGGUUAG'back_transcribe(my_string)# 'GCTGTTATGGGTCGTTGGAAGGGTGGTCGTGCTGCTGGTTAG'translate(my_string)# 'AVMGRWKGGRAAG*'

实际工作中，序列往往伴随着大量注释信息（如基因位置、功能描述、参考文献等）。Biopython 使用 SeqRecord 对象将序列和注释整合在一起。

SeqRecord 定义在 Bio.SeqRecord 模块中，主要包含以下属性：

• .seq：序列本身，通常是 Seq 对象。

• .id：主标识符（如 accession number）。

• .name：常用名称（类似 GenBank 的 LOCUS 名）。

• .description：人类可读的描述。

• .letter_annotations：每个字母的注释字典（如质量分数、二级结构）。

• .annotations：关于整个序列的注释字典（如来源、关键词）。

• .features：SeqFeature 对象列表，描述序列上的特征（如基因、CDS）。

• .dbxrefs：数据库交叉引用列表（如“Project:58037”）。

from Bio.SeqRecord import SeqRecordhelp(SeqRecord)

最简单的方式是直接从 Seq 对象创建，并指定必要的信息：

from Bio.Seq import Seqfrom Bio.SeqRecord import SeqRecordsimple_seq = Seq("GATC")simple_seq_r = SeqRecord(simple_seq)simple_seq_r.id    # '<unknown id>'record = SeqRecord(simple_seq, id="AC12345", name="AC12345",                   description="A random sequence")record.id    # 'AC12345'

simple_seq_r.id = "AC12345"simple_seq_r.description = "Made up sequence I wish I could write a paper about"

通常不会手动创建 SeqRecord，而是用 Bio.SeqIO 从文件读取。

1. 从 FASTA 文件读取

# 文件下载链接 https://raw.githubusercontent.com/biopython/biopython/master/Tests/GenBank/NC_005816.fnafrom Bio import SeqIOrecord = SeqIO.read("NC_005816.fna", "fasta")print(record.id)          # gi|45478711|ref|NC_005816.1|print(record.description) # 整个标题行print(record.seq)         # 序列对象print(record.annotations) # 空字典（FASTA 没有注释）print(record.features)    # 空列表

2. 从 GenBank 文件读取

GenBank 文件包含丰富的注释，读取后 annotations 和 features 都会被填充：

# 文件下载链接 https://raw.githubusercontent.com/biopython/biopython/master/Tests/GenBank/NC_005816.gbfrom Bio import SeqIOrecord = 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'])

特征 (feature) 是序列上的一个区域，通常带有类型和位置信息。SeqFeature 的主要属性：

• .type：特征类型（字符串），如“gene”、“CDS”。

• .location：位置对象，描述区域。

• .qualifiers：附加信息的字典（如基因名、产物），值总是字符串列表。

• .strand：链的方向（+1、-1、0 或 None）。

位置 (location) 可以简单（如 5..20），也可以复杂（如 join 或模糊位置）。Biopython 使用 SimpleLocation 和 CompoundLocation 表示，位置本身可以是精确或模糊的（如 <13、>13、(1.5) 等）。

例如，创建一个带有模糊起止点的位置：

from Bio import SeqFeaturestart_pos = SeqFeature.AfterPosition(5)end_pos = SeqFeature.BetweenPosition(9, left=8, right=9)my_location = SeqFeature.SimpleLocation(start_pos, end_pos)print(my_location)# [>5:(8^9)]

SeqFeature 对象本身不包含序列，但可以用 .extract() 方法从父序列中提取：

from Bio.Seq import Seqfrom Bio.SeqFeature import SeqFeature, SimpleLocationseq = Seq("ACCGAGACGGCAAAGGCTAGCATAGGTATGAGACTTCCTTCCTGCCAGTGCTGAGGAACTGGGAGCCTAC")feature = SeqFeature(SimpleLocation(5, 18, strand=-1), type="gene")feature_seq = feature.extract(seq)feature_seq# Seq('AGCCTTTGCCGTC')# 对父序列进行切片以提取 5:18, 然后取反向互补# 以下代码也是一样的效果feature_seq = seq[feature.location.start : feature.location.end].reverse_complement()print(feature_seq)# AGCCTTTGCCGTC

直接比较两个 SeqRecord 对象会抛出 NotImplementedError，因为比较逻辑太复杂。建议只比较你关心的属性：

record1 = SeqRecord(Seq("ACGT"), id="test")record2 = SeqRecord(Seq("ACGT"), id="test")print(record1.id == record2.id)   # Trueprint(record1.seq == record2.seq) # True

Bio.SeqFeature.Reference 类可以存储参考文献信息，例如标题、作者、期刊、PubMed ID 等。这些 Reference 对象通常放在 SeqRecord.annotations["references"] 列表中。

SeqRecord 有一个 .format() 方法，可以将记录输出为指定格式的字符串（如 FASTA、GenBank）：

print(record.format("fasta"))

注意：某些格式需要多个记录才能写入（如多序列比对格式），此时 .format() 会失败。

对 SeqRecord 切片会得到一个新的 SeqRecord，包含对应区域的序列。特征和每个字母的注释也会被相应调整，但大部分顶级注释（如 annotations 和 dbxrefs）会被丢弃，以避免错误传播。你可以手动恢复需要的注释。

from Bio import SeqIOrecord = SeqIO.read("NC_005816.gb", "genbank")record# SeqRecord(seq=Seq('TGTAACGAACGGTGCAAT...CTG'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=['Project:58037'])len(record)# 9609len(record.features)# 41sub_record = record[4300:4800]sub_record# SeqRecord(seq=Seq('ATAAATAGATTATTCCAAA...TTA'), id='NC_005816.1', name='NC_005816', description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])len(sub_record)# 500len(sub_record.features)# 2

两个 SeqRecord 可以用 + 拼接，生成一个新的记录。每个字母的注释会被拼接，特征会被调整到新坐标，公共注释（如 id、name）会被保留，但 annotations 和 dbxrefs 仅当两个记录相同时才保留。

例如，从 FASTQ 文件中拼接两个片段以纠正一个碱基：

record = next(SeqIO.parse("example.fastq", "fastq"))left = record[:20]right = record[21:]edited = left + right   # 跳过第21个碱基

SeqRecord 也有 .reverse_complement() 方法，它会返回一个反向互补的新记录。序列本身被反向互补，特征的位置和链被调整，每个字母注释被反转。但默认情况下，id、name、description、annotations 和 dbxrefs 不会被复制，你需要显式指定是否保留或提供新值。

from Bio import SeqIOrec = 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 13rc = rec.reverse_complement(id="TESTING")print(rc.id, len(rc), len(rc.features), len(rc.dbxrefs), len(rc.annotations))# TESTING 9609 41 0 0

本篇详细介绍了 Biopython 的两个核心对象：

• Seq：轻量级的序列容器，提供了丰富的生物学方法（互补、转录、翻译等）。

• SeqRecord：扩展了 Seq，集成了序列的所有注释信息（标识符、特征、参考文献等），是序列输入输出（Bio.SeqIO）的基础。