biopython-sequence-analysis

Biopython：序列分析工具包

概述

Biopython 提供了一套全面的模块，用于以序列为中心的生物信息学工作：读取和写入各种主要的生物学文件格式（FASTA、FASTQ、GenBank、GFF），以编程方式查询 NCBI 数据库，运行 BLAST 搜索并解析结果，对序列进行两两比对或多序列比对，以及构建和可视化系统发育树。本技能专注于分析工作流——从检索 NCBI 数据，到序列比对，再到系统发育推断。

如需了解 PCR 引物设计、限制性内切酶消化、克隆模拟、蛋白质结构分析（Bio.PDB）以及分子量/Tm 计算，请参阅 biopython-molecular-biology。

适用场景

• 从 NCBI Nucleotide/Protein 下载一个基因家族，对序列进行比对并构建系统发育树
• 解析 GenBank 或 GFF3 注释文件，并提取一组特征的 CDS 序列
• 对 NCBI nt 或 nr 运行 BLAST 搜索，筛选显著命中项并获取其完整序列
• 使用 BLOSUM62/PAM250 矩阵计算两两序列一致性或为比对结果评分
• 使用 SeqIO.index() 为大型多序列 FASTA 或 FASTQ 文件创建索引，以便随机访问检索，而无需将所有序列加载到 RAM
• 通过一次调用在不同序列格式间转换（FASTA ↔ GenBank ↔ FASTQ ↔ PHYLIP）
• 以编程方式遍历、定根、修剪和注释 Newick 或 Nexus 系统发育树
• 处理 SAM/BAM/CRAM 比对文件和比对上的 reads 时，请改用 pysam
• 需要生态多样性指标（UniFrac、beta 多样性）或 PCoA 等排序方法时，请改用 scikit-bio
• 无需编写 Entrez 工作流、只需快速查询基因或通过一行命令查询 NCBI 时，请改用 gget

先决条件

• Python 包：biopython、numpy、matplotlib
• 可选工具：本地安装 MUSCLE 或 ClustalW（用于 Bio.Align.Applications 包装器）
• NCBI 访问：在调用任何 E-utilities 前设置 Entrez.email；前往 https://www.ncbi.nlm.nih.gov/account/ 获取免费 API 密钥，可将请求速率提升至 10 次/秒（默认速率为 3 次/秒）
• 本地 BLAST：需单独安装 BLAST+（离线搜索可运行 conda install -c bioconda blast）

安装前请先检查：当前环境中可能已经提供该工具（例如，在 pixi / conda 环境中）。先运行 command -v python；如果返回路径，请跳过下方的安装命令。在 pixi 项目中运行时，请使用 pixi run python，而不是直接运行 python。

pip install biopython numpy matplotlib
conda install -c bioconda blast  # optional, for local BLAST

快速入门

from Bio import SeqIO, Entrez
from Bio.SeqUtils import gc_fraction

# Fetch a GenBank record and display basic stats
Entrez.email = "[email protected]"
handle = Entrez.efetch(db="nucleotide", id="NM_007294", rettype="gb", retmode="text")
record = SeqIO.read(handle, "genbank")
handle.close()

print(f"ID: {record.id}")
print(f"Length: {len(record.seq)} bp")
print(f"GC content: {gc_fraction(record.seq)*100:.1f}%")
print(f"Features: {len(record.features)}")
print(f"First feature: {record.features[0].type} at {record.features[0].location}")
# ID: NM_007294.4
# Length: 7207 bp
# GC content: 47.3%
# Features: 9

核心 API

模块 1：SeqIO — 文件解析和格式转换

Bio.SeqIO 可读取和写入各种主要的序列格式（FASTA、FASTQ、GenBank、EMBL、PHYLIP、Nexus）。SeqIO.parse() 返回由 SeqRecord 对象组成的迭代器；SeqIO.index() 则会为大型文件构建磁盘或内存字典。

from Bio import SeqIO

# Parse FASTA and FASTQ
fasta_records = list(SeqIO.parse("sequences.fasta", "fasta"))
print(f"FASTA: {len(fasta_records)} sequences")

# FASTQ: access quality scores
for rec in SeqIO.parse("reads.fastq", "fastq"):
    quals = rec.letter_annotations["phred_quality"]
    avg_q = sum(quals) / len(quals)
    print(f"  {rec.id}: {len(rec.seq)} bp, mean Q={avg_q:.1f}")
    break  # show one example

# Parse GenBank with feature access
for rec in SeqIO.parse("chromosome.gb", "genbank"):
    cdss = [f for f in rec.features if f.type == "CDS"]
    print(f"{rec.id}: {len(cdss)} CDS features")
    for cds in cdss[:3]:
        gene = cds.qualifiers.get("gene", ["unknown"])[0]
        print(f"  {gene}: {cds.location}")

from Bio import SeqIO

# SeqIO.index() for random access without loading all records
# Useful for large reference FASTA files (genomes, nr database subsets)
idx = SeqIO.index("large_genome.fasta", "fasta")
print(f"Index contains {len(idx)} sequences")

# Retrieve specific sequences by ID in O(1)
target = idx["chr1"]
print(f"chr1: {len(target.seq):,} bp")
region = target.seq[1_000_000:1_001_000]
print(f"Region [1M-1M+1kb]: {region[:60]}...")

# Format conversion: GenBank → FASTA in one call
n = SeqIO.convert("annotation.gb", "genbank", "sequences.fasta", "fasta")
print(f"Converted {n} records to FASTA")

模块 2：Seq 和 SeqRecord — 序列对象与特征注释

Bio.Seq 以支持原地操作的方式表示生物学序列。Bio.SeqRecord 将 Seq 封装起来，并附带 ID、描述以及特征注释字典。特征使用带链向信息（+1、-1）的 FeatureLocation。

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio.SeqUtils import gc_fraction, MeltingTemp

dna = Seq("ATGAAACCCGGGTTTTAA")
print(f"GC content: {gc_fraction(dna)*100:.1f}%")
print(f"Reverse complement: {dna.reverse_complement()}")
print(f"Transcript: {dna.transcribe()}")
print(f"Translation: {dna.translate()}")
print(f"Translation (to stop): {dna.translate(to_stop=True)}")

# Tm calculation for a short primer
primer = Seq("ATGAAACCCGGG")
tm = MeltingTemp.Tm_Wallace(primer)
print(f"Tm (Wallace): {tm:.1f}°C")

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation

# Build a SeqRecord with annotated features
gene_seq = Seq("ATGAAACCCGGGTTTTAAATCGATCG" * 10)
record = SeqRecord(gene_seq, id="MY_GENE_001", description="synthetic example gene")

# Annotate a CDS feature
cds = SeqFeature(
    FeatureLocation(0, 18, strand=+1),
    type="CDS",
    qualifiers={"gene": ["myGene"], "product": ["hypothetical protein"]}
)
record.features.append(cds)

# Extract and translate the feature
cds_seq = cds.location.extract(record.seq)
protein = cds_seq.translate(to_stop=True)
print(f"CDS: {cds_seq}")
print(f"Protein: {protein}")

# Slice record preserving feature annotations
sub = record[0:60]
print(f"Subrecord: {len(sub.seq)} bp, {len(sub.features)} features")

模块 3：Entrez — 以编程方式访问 NCBI 数据库

Bio.Entrez 封装了 NCBI E-utilities API：esearch 查找符合查询条件的记录，efetch 获取完整记录，elink 查找不同数据库间的相关记录，esummary 返回文档摘要。务必设置 Entrez.email，并遵守每秒 3 次请求的速率限制（使用 API 密钥时为每秒 10 次）。

from Bio import Entrez, SeqIO
import time

Entrez.email = "[email protected]"
# Entrez.api_key = "YOUR_API_KEY"  # for 10 req/s

# esearch: find IDs matching a query
handle = Entrez.esearch(db="nucleotide", term="BRCA1[Gene] AND Homo sapiens[Organism] AND mRNA[Filter]", retmax=10)
search_results = Entrez.read(handle)
handle.close()

print(f"Total hits: {search_results['Count']}")
ids = search_results["IdList"]
print(f"Retrieved IDs: {ids}")

# efetch: download full GenBank records
for acc_id in ids[:3]:
    handle = Entrez.efetch(db="nucleotide", id=acc_id, rettype="gb", retmode="text")
    record = SeqIO.read(handle, "genbank")
    handle.close()
    print(f"  {record.id}: {len(record.seq)} bp — {record.description[:60]}")
    time.sleep(0.4)  # stay within rate limit

from Bio import Entrez
import time

Entrez.email = "[email protected]"

# elink: cross-database links (e.g., PubMed article → related nucleotide sequences)
handle = Entrez.elink(dbfrom="pubmed", db="nucleotide", id="29087512")
link_results = Entrez.read(handle)
handle.close()

linked_ids = []
for linkset in link_results:
    for db_links in linkset.get("LinkSetDb", []):
        if db_links["DbTo"] == "nucleotide":
            linked_ids = [lnk["Id"] for lnk in db_links["Link"]]
            break

print(f"Nucleotide sequences linked to PubMed 29087512: {len(linked_ids)}")
print(f"First IDs: {linked_ids[:5]}")

# esummary: lightweight metadata without downloading full records
if linked_ids:
    handle = Entrez.esummary(db="nucleotide", id=",".join(linked_ids[:5]))
    summaries = Entrez.read(handle)
    handle.close()
    for doc in summaries:
        print(f"  {doc['AccessionVersion']}: {doc['Title'][:70]}")

模块 4：BLAST — 远程和本地序列相似性搜索

Bio.Blast.NCBIWWW.qblast() 会向 NCBI BLAST 服务提交查询并返回 XML 句柄。Bio.Blast.NCBIXML.parse()（单个查询时可使用 .read()）会生成包含 alignments 和 hsps 的 Blast 对象。进行大规模搜索时，请通过 subprocess 使用本地 BLAST+。

from Bio.Blast import NCBIWWW, NCBIXML
from Bio.Seq import Seq

# Remote BLASTP against Swiss-Prot (small, reviewed database)
query = Seq("MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSY")
result_handle = NCBIWWW.qblast("blastp", "swissprot", str(query), hitlist_size=10)

# Parse results
blast_record = NCBIXML.read(result_handle)
print(f"Query: {blast_record.query}")
print(f"Database: {blast_record.database}")
print(f"Hits: {len(blast_record.alignments)}")

E_VALUE_THRESH = 1e-5
for alignment in blast_record.alignments[:5]:
    for hsp in alignment.hsps:
        if hsp.expect < E_VALUE_THRESH:
            identity_pct = hsp.identities / hsp.align_length * 100
            print(f"\n  Hit: {alignment.title[:70]}")
            print(f"  E-value: {hsp.expect:.2e}, Identity: {identity_pct:.1f}%, Score: {hsp.score}")
            print(f"  Query:  {hsp.query[:60]}")
            print(f"  Match:  {hsp.match[:60]}")
            print(f"  Sbjct:  {hsp.sbjct[:60]}")

import subprocess
from Bio.Blast import NCBIXML
from Bio import SeqIO

# Local BLAST+ via subprocess (faster for large batches)
# Requires: makeblastdb and blastp installed (conda install -c bioconda blast)

# Step 1: Build a local database from a FASTA file
subprocess.run(
    ["makeblastdb", "-in", "ref_proteins.fasta", "-dbtype", "prot", "-out", "ref_db"],
    check=True
)

# Step 2: Run blastp against local database
result = subprocess.run(
    ["blastp", "-query", "query.fasta", "-db", "ref_db",
     "-outfmt", "5",           # XML output for NCBIXML parsing
     "-evalue", "1e-5",
     "-num_threads", "4",
     "-out", "blast_results.xml"],
    check=True
)

# Step 3: Parse XML results
with open("blast_results.xml") as fh:
    for blast_rec in NCBIXML.parse(fh):
        print(f"Query: {blast_rec.query_id}")
        for aln in blast_rec.alignments[:3]:
            hsp = aln.hsps[0]
            print(f"  Hit: {aln.hit_id}, E={hsp.expect:.2e}, Id={hsp.identities}/{hsp.align_length}")

模块 5：Phylo — 树解析、操作与可视化

Bio.Phylo 可读取 Newick、Nexus、PhyloXML 和 NeXML 格式。树由带有嵌套子节点的 Clade 对象表示。主要方法：find_clades()（遍历节点的生成器）、common_ancestor()、distance()、root_with_outgroup() 和 prune()。可视化使用 matplotlib。

from Bio import Phylo
import io

# Parse a Newick tree from a string
newick = "((Homo_sapiens:0.01, Pan_troglodytes:0.012):0.08, (Mus_musculus:0.15, Rattus_norvegicus:0.14):0.12, Drosophila_melanogaster:0.85);"
tree = Phylo.read(io.StringIO(newick), "newick")

print(f"Number of terminals: {tree.count_terminals()}")
print(f"Total branch length: {tree.total_branch_length():.3f}")
print(f"Terminals: {[t.name for t in tree.get_terminals()]}")

# Root with outgroup and compute distances
tree.root_with_outgroup("Drosophila_melanogaster")
human = tree.find_any("Homo_sapiens")
chimp = tree.find_any("Pan_troglodytes")
mouse = tree.find_any("Mus_musculus")
print(f"\nDistance Human–Chimp: {tree.distance(human, chimp):.4f}")
print(f"Distance Human–Mouse: {tree.distance(human, mouse):.4f}")

# Traverse all internal nodes
for clade in tree.find_clades(order="level"):
    if not clade.is_terminal():
        children = [c.name or "internal" for c in clade.clades]
        print(f"  Internal node → {children}")

from Bio import Phylo
import matplotlib.pyplot as plt
import io

newick = "((Homo_sapiens:0.01, Pan_troglodytes:0.012):0.08, (Mus_musculus:0.15, Rattus_norvegicus:0.14):0.12, Drosophila_melanogaster:0.85);"
tree = Phylo.read(io.StringIO(newick), "newick")
tree.root_with_outgroup("Drosophila_melanogaster")
tree.ladderize()   # sort clades by size for a clean layout

# Annotate bootstrap support (if present in node labels)
for clade in tree.find_clades():
    if clade.confidence is not None:
        clade.name = f"{clade.confidence:.0f}"

fig, ax = plt.subplots(figsize=(8, 5))
Phylo.draw(tree, axes=ax, do_show=False, show_confidence=True)
ax.set_title("Primate + Outgroup Phylogeny")
plt.tight_layout()
plt.savefig("phylogeny.png", dpi=150, bbox_inches="tight")
print("Saved phylogeny.png")

模块 6：PairwiseAligner — 两两序列比对

Bio.Align.PairwiseAligner 取代了旧版 pairwise2 模块（自 Biopython 1.80 起弃用）。它支持局部和全局比对、缺口开放／延伸罚分，以及 Bio.Align.substitution_matrices 中的任意替换矩阵。

from Bio.Align import PairwiseAligner, substitution_matrices

# Global protein alignment with BLOSUM62
aligner = PairwiseAligner()
aligner.mode = "global"
aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
aligner.open_gap_score = -11
aligner.extend_gap_score = -1

seq1 = "MTEYKLVVVGAGGVGKSALTIQLIQNHFVDEYDPTIEDSY"
seq2 = "MTEYKLVVVGAVGVGKSALTIQLIQNHFVDEYDPTIEDSY"  # G12V variant

alignments = aligner.align(seq1, seq2)
best = alignments[0]
print(f"Score: {best.score:.1f}")
print(f"Number of alignments: {aligner.score(seq1, seq2):.0f} (score only, faster)")
print(best)   # formatted alignment string

identity = sum(a == b for a, b in zip(*best) if a != "-") / best.shape[1]
print(f"Identity: {identity*100:.1f}%")

from Bio.Align import PairwiseAligner, substitution_matrices

# Local DNA alignment
aligner = PairwiseAligner()
aligner.mode = "local"
aligner.match_score = 2
aligner.mismatch_score = -1
aligner.open_gap_score = -5
aligner.extend_gap_score = -0.5

query = "ACGTACGTACGT"
subject = "TTTTACGTACGTACGTTTTT"
alignments = list(aligner.align(query, subject))
print(f"Local alignments found: {len(alignments)}")
best = alignments[0]
print(f"Score: {best.score}")
print(f"Aligned region in subject: [{best.coordinates[1][0]}:{best.coordinates[1][-1]}]")
print(best)

# Batch pairwise identity matrix
seqs = ["ACGTACGT", "ACGTATGT", "TTGTACGT", "ACGTACGG"]
n = len(seqs)
aligner2 = PairwiseAligner(mode="global", match_score=1, mismatch_score=-1)
print("\nPairwise identity matrix:")
for i in range(n):
    row = []
    for j in range(n):
        score = aligner2.score(seqs[i], seqs[j])
        max_len = max(len(seqs[i]), len(seqs[j]))
        row.append(f"{score/max_len:.2f}")
    print("  " + "  ".join(row))

核心概念

SeqRecord 和特征坐标

SeqRecord 用于存储序列元数据和 SeqFeature 对象列表。特征使用从零开始、左闭右开的 FeatureLocation(start, end, strand)，与 Python 切片使用相同的约定。使用 feature.location.extract(record.seq) 获取链方向正确的子序列。复合位置（例如剪接后的外显子）使用 CompoundLocation。

from Bio import SeqIO

# Extract all CDS protein translations from a GenBank record
for rec in SeqIO.parse("gene.gb", "genbank"):
    for feat in rec.features:
        if feat.type == "CDS":
            gene = feat.qualifiers.get("gene", ["?"])[0]
            product = feat.qualifiers.get("product", ["unknown"])[0]
            cds_nt = feat.location.extract(rec.seq)
            protein = cds_nt.translate(to_stop=True)
            print(f"{gene} ({product}): {len(protein)} aa — {protein[:10]}...")

Phylo Clade 对象

Clade 既是一个节点，也是以该节点为根的子树。末端 clade（叶节点）设有 .name；内部 clade 可能设有 .confidence（自举值）和 .branch_length。使用 tree.find_clades() 并设置 terminal=True/False 进行筛选。根节点是 tree.root（同样是一个 Clade）。Phylo.read() 返回一个 Tree 包装对象；tree.root 返回根 Clade。

常见工作流

工作流 1：基因家族系统发育分析 — NCBI 下载 → MUSCLE 比对 → 构建树

目标：下载脊椎动物中的所有 BRCA1 直系同源基因，使用 MUSCLE 进行比对，并构建邻接法树。

import subprocess
import time
from Bio import Entrez, SeqIO, Phylo, AlignIO
from Bio.Align import MultipleSeqAlignment
import matplotlib.pyplot as plt

Entrez.email = "[email protected]"

# Step 1: Search NCBI Protein for BRCA1 orthologs
handle = Entrez.esearch(
    db="protein",
    term="BRCA1[Gene] AND Vertebrata[Organism] AND RefSeq[Filter]",
    retmax=20
)
results = Entrez.read(handle); handle.close()
ids = results["IdList"]
print(f"Found {results['Count']} sequences, using {len(ids)}")

# Step 2: Fetch sequences in FASTA format
handle = Entrez.efetch(db="protein", id=",".join(ids), rettype="fasta", retmode="text")
with open("brca1_orthologs.fasta", "w") as fh:
    fh.write(handle.read())
handle.close()
print(f"Saved {len(ids)} sequences to brca1_orthologs.fasta")
time.sleep(1)

# Step 3: Align with MUSCLE (must be installed: conda install -c bioconda muscle)
subprocess.run(
    ["muscle", "-align", "brca1_orthologs.fasta", "-output", "brca1_aligned.fasta"],
    check=True
)
alignment = AlignIO.read("brca1_aligned.fasta", "fasta")
print(f"Alignment: {len(alignment)} sequences × {alignment.get_alignment_length()} columns")

# Step 4: Build a neighbor-joining tree using Bio.Phylo + distance matrix
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor

calculator = DistanceCalculator("identity")
dm = calculator.get_distance(alignment)
constructor = DistanceTreeConstructor(calculator, method="nj")
tree = constructor.build_tree(alignment)

# Step 5: Root and save
tree.root_at_midpoint()
tree.ladderize()
Phylo.write(tree, "brca1_nj.nwk", "newick")
print("Saved NJ tree to brca1_nj.nwk")

# Step 6: Visualize
fig, ax = plt.subplots(figsize=(10, 8))
Phylo.draw(tree, axes=ax, do_show=False)
ax.set_title("BRCA1 Ortholog NJ Tree")
plt.tight_layout()
plt.savefig("brca1_tree.png", dpi=150, bbox_inches="tight")
print("Saved brca1_tree.png")

工作流 2：BLAST 流程 — FASTA 查询 → 远程 BLAST → 获取排名靠前的匹配项

目标：读取 FASTA 查询序列，在 NCBI nr 数据库中进行 BLAST，筛选显著匹配项，获取排名靠前的匹配项的完整 GenBank 记录，并写入汇总 CSV 文件。

import csv
import time
from Bio import SeqIO, Entrez
from Bio.Blast import NCBIWWW, NCBIXML

Entrez.email = "[email protected]"

# Step 1: Load query sequence
query_record = SeqIO.read("query.fasta", "fasta")
print(f"Query: {query_record.id} ({len(query_record.seq)} aa)")

# Step 2: Remote BLAST against nr protein database
print("Submitting BLAST job (may take 1-5 minutes)...")
result_handle = NCBIWWW.qblast(
    "blastp", "nr", str(query_record.seq),
    hitlist_size=20,
    expect=1e-5,
    word_size=6,
    matrix_name="BLOSUM62"
)

# Step 3: Parse results and filter
E_THRESH = 1e-5
MIN_IDENTITY = 50.0
top_hits = []

blast_record = NCBIXML.read(result_handle)
for alignment in blast_record.alignments:
    for hsp in alignment.hsps:
        if hsp.expect > E_THRESH:
            continue
        identity_pct = hsp.identities / hsp.align_length * 100
        if identity_pct < MIN_IDENTITY:
            continue
        accession = alignment.accession
        top_hits.append({
            "accession": accession,
            "title": alignment.title[:80],
            "length": alignment.length,
            "score": hsp.score,
            "evalue": hsp.expect,
            "identity_pct": round(identity_pct, 1),
            "coverage": round(hsp.align_length / blast_record.query_length * 100, 1),
        })

print(f"Filtered to {len(top_hits)} significant hits")

# Step 4: Fetch full GenBank records for top 5 hits
accessions = [h["accession"] for h in top_hits[:5]]
time.sleep(1)
handle = Entrez.efetch(db="protein", id=",".join(accessions), rettype="gb", retmode="text")
gb_records = list(SeqIO.parse(handle, "genbank"))
handle.close()
SeqIO.write(gb_records, "top_blast_hits.gb", "genbank")
print(f"Saved {len(gb_records)} full GenBank records to top_blast_hits.gb")

# Step 5: Write CSV summary
with open("blast_summary.csv", "w", newline="") as fh:
    writer = csv.DictWriter(fh, fieldnames=top_hits[0].keys())
    writer.writeheader()
    writer.writerows(top_hits)
print(f"Saved blast_summary.csv ({len(top_hits)} hits)")

工作流 3：FASTQ 质量筛选与格式转换流程

目标：读取 FASTQ 文件，根据平均质量和长度筛选 reads，并输出供下游比对使用的 FASTA 文件。

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

input_fastq = "raw_reads.fastq"
output_fasta = "filtered_reads.fasta"
MIN_QUAL = 20     # mean Phred quality
MIN_LEN = 100     # minimum read length
MAX_LEN = 300     # maximum read length

def passes_qc(record, min_q=MIN_QUAL, min_l=MIN_LEN, max_l=MAX_LEN):
    quals = record.letter_annotations["phred_quality"]
    mean_q = sum(quals) / len(quals)
    length = len(record.seq)
    return mean_q >= min_q and min_l <= length <= max_l

kept = 0
total = 0
with open(output_fasta, "w") as out_fh:
    for record in SeqIO.parse(input_fastq, "fastq"):
        total += 1
        if passes_qc(record):
            # Convert FASTQ → FASTA (drops quality scores)
            SeqIO.write(record, out_fh, "fasta")
            kept += 1

print(f"Total reads: {total:,}")
print(f"Passed QC: {kept:,} ({kept/total*100:.1f}%)")
print(f"Saved to: {output_fasta}")

关键参数

参数	模块／函数	默认值	范围／选项	作用
retmax	Entrez.esearch()	20	1–10000	每次搜索返回的最大 ID 数；超过 10000 时使用历史服务器
hitlist_size	NCBIWWW.qblast()	50	1–5000	返回的最大 BLAST 比对数
expect	NCBIWWW.qblast()	10.0	1e-100–1000	BLAST 匹配项报告的 E 值阈值
matrix_name	NCBIWWW.qblast()	"BLOSUM62"	"BLOSUM45"、"BLOSUM80"、"PAM250"	替换矩阵；常规用途使用 BLOSUM62
mode	PairwiseAligner	"global"	"global"、"local"	Needleman-Wunsch（全局）或 Smith-Waterman（局部）
open_gap_score	PairwiseAligner	-1	负浮点数	开启缺口的罚分；增大绝对值可提高罚分
extend_gap_score	PairwiseAligner	0	负浮点数	每个残基的缺口延伸罚分
method	DistanceTreeConstructor	"nj"	"nj"、"upgma"	NJ（无根）或 UPGMA（有根，假定分子钟）
format	Phylo.read/write()	必需	"newick"、"nexus"、"phyloxml"	树文件格式
rettype	Entrez.efetch()	必需	"fasta"、"gb"、"xml"	返回的记录格式；需搭配相应的 SeqIO 格式字符串

最佳实践

1. 始终设置 Entrez.email 并遵守速率限制。NCBI 会封禁每秒请求超过 3 次且未使用密钥的 IP。在循环中的请求之间添加 time.sleep(0.4)，或使用 Entrez.api_key 将速率提高到每秒 10 次。

   from Bio import Entrez
   import time
   Entrez.email = "[email protected]"
   Entrez.api_key = "YOUR_API_KEY"    # from https://www.ncbi.nlm.nih.gov/account/
   # In any batch loop: time.sleep(0.11)  # ~9 req/s with key
   

2. 处理大型 FASTA 文件时使用 SeqIO.index()，而不是 list(SeqIO.parse(...))。将所有记录加载到列表中会占用 O(N) 内存；index() 只读取字节偏移量，并按需获取记录。

   # Bad for large files: loads everything into RAM
   # records = {r.id: r for r in SeqIO.parse("genome.fasta", "fasta")}
   
   # Good: on-disk index, O(1) access
   from Bio import SeqIO
   idx = SeqIO.index("genome.fasta", "fasta")
   rec = idx["chr22"]          # fetches only this record
   

3. 使用 PairwiseAligner，不要使用旧版 pairwise2。自 Biopython 1.80 起，Bio.pairwise2 已弃用，并将被移除。PairwiseAligner 速度更快，支持直接使用替换矩阵，并返回带坐标数组的 Alignment 对象。

4. 读取后立即关闭 Entrez 句柄。句柄是 HTTP 连接；让它们保持打开状态可能导致超时和资源耗尽。

   handle = Entrez.efetch(db="nucleotide", id="NM_007294", rettype="gb", retmode="text")
   record = SeqIO.read(handle, "genbank")
   handle.close()    # do not skip this
   

5. 批量搜索时优先使用本地 BLAST。远程 NCBIWWW.qblast() 适用于临时查询，但在 NCBI 服务器上可能排队数分钟。筛查超过 100 条序列时，使用 makeblastdb 构建本地 BLAST+ 数据库，并通过 subprocess 调用 blastp/blastn，使用 -outfmt 5（XML）以便用 NCBIXML 解析。

6. 测量距离前先为树定根。Phylo.distance() 计算两个节点之间路径上的分支长度总和。在无根树上，路径仍然是唯一的，但对于可视化和分支断言，使用中点定根或外群定根在生物学上更合理。

7. 构建 SeqRecord 集合前，先去除比对中的间隙。BLAST 和比对结果可能包含间隙字符（-）。对含间隙序列执行 .translate() 等 Seq 操作会报错；请使用 seq.replace("-", "") 去除间隙，或使用 ungap()。

常用示例

示例：使用 History Server 批量获取 Entrez 记录

适用场景：从 NCBI 下载超过 500 条记录，可避免 URL 长度限制，并将搜索结果保留在服务器端。

from Bio import Entrez, SeqIO
import time

Entrez.email = "[email protected]"

# Search and store results on NCBI history server
handle = Entrez.esearch(db="nucleotide", term="16S rRNA[Gene] AND Bacteria[Organism]",
                        retmax=500, usehistory="y")
results = Entrez.read(handle); handle.close()
webenv = results["WebEnv"]
query_key = results["QueryKey"]
count = int(results["Count"])
print(f"Total: {count} sequences")

# Fetch in batches of 200
batch_size = 200
all_records = []
for start in range(0, min(count, 1000), batch_size):
    handle = Entrez.efetch(
        db="nucleotide", rettype="fasta", retmode="text",
        retstart=start, retmax=batch_size,
        webenv=webenv, query_key=query_key
    )
    batch = list(SeqIO.parse(handle, "fasta"))
    handle.close()
    all_records.extend(batch)
    print(f"  Fetched {len(all_records)}/{min(count, 1000)}")
    time.sleep(0.5)

SeqIO.write(all_records, "16S_sequences.fasta", "fasta")
print(f"Saved {len(all_records)} sequences")

示例：解析带序列特征的 GFF3

适用场景：从与参考 FASTA 配对的 GFF3 注释中提取基因/CDS 序列。

from Bio import SeqIO

# Load genome FASTA into indexed dict
genome = SeqIO.to_dict(SeqIO.parse("genome.fasta", "fasta"))

# Parse GFF3 manually (Biopython does not have a native GFF3 parser;
# use the gffutils package for complex queries, or parse directly for simple cases)
genes = []
with open("annotation.gff3") as fh:
    for line in fh:
        if line.startswith("#") or not line.strip():
            continue
        fields = line.strip().split("\t")
        if len(fields) < 9 or fields[2] != "CDS":
            continue
        chrom, _, feat_type, start, end, _, strand, _, attrs = fields
        start, end = int(start) - 1, int(end)   # GFF3 is 1-based
        if chrom not in genome:
            continue
        seq = genome[chrom].seq[start:end]
        if strand == "-":
            seq = seq.reverse_complement()
        attr_dict = dict(kv.split("=") for kv in attrs.strip().split(";") if "=" in kv)
        gene_id = attr_dict.get("gene_id", attr_dict.get("ID", "unknown"))
        genes.append((gene_id, seq, strand))

print(f"Extracted {len(genes)} CDS features")
for gid, seq, strand in genes[:3]:
    print(f"  {gid} ({strand}): {len(seq)} bp — protein: {seq.translate(to_stop=True)[:8]}...")

示例：根据多序列比对计算成对一致性矩阵

适用场景：快速评估比对文件（FASTA、PHYLIP、Clustal）中的序列多样性，并识别离群序列。

from Bio import AlignIO
import numpy as np

alignment = AlignIO.read("aligned_sequences.fasta", "fasta")
n = len(alignment)
names = [rec.id for rec in alignment]
length = alignment.get_alignment_length()

# Build identity matrix
identity_matrix = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        if i == j:
            identity_matrix[i, j] = 1.0
            continue
        matches = sum(
            a == b and a != "-"
            for a, b in zip(str(alignment[i].seq), str(alignment[j].seq))
        )
        aligned_cols = sum(a != "-" and b != "-"
                           for a, b in zip(str(alignment[i].seq), str(alignment[j].seq)))
        identity_matrix[i, j] = matches / aligned_cols if aligned_cols else 0.0

print("Pairwise identity matrix:")
print(f"{'':20s} " + "  ".join(f"{n[:8]:>8s}" for n in names))
for i, name in enumerate(names):
    row = "  ".join(f"{identity_matrix[i,j]*100:8.1f}" for j in range(n))
    print(f"{name[:20]:20s} {row}")

# Find most divergent pair
min_id = np.min(identity_matrix[identity_matrix > 0])
idx = np.unravel_index(np.argmin(np.where(identity_matrix > 0, identity_matrix, 1)), identity_matrix.shape)
print(f"\nMost divergent pair: {names[idx[0]]} vs {names[idx[1]]} ({min_id*100:.1f}%)")

故障排除

问题	原因	解决方案
urllib.error.HTTPError: 429 Too Many Requests	超出 NCBI 速率限制（3 次请求/秒）	在调用之间添加 time.sleep(0.4)，或在 https://www.ncbi.nlm.nih.gov/account/ 注册免费 API 密钥，将速率提升至 10 次请求/秒
RuntimeError: Too many requests were made without pausing	NCBI 检测到连续快速请求	设置 Entrez.email（必需），降低循环频率，将 ID 合并为逗号分隔的字符串，并在一次 efetch 调用中批量处理
Bio.Application.ApplicationError: blastall not found	旧版 blastall 已被 BLAST+ 取代	将 NcbiblastpCommandline 替换为 subprocess.run(["blastp", ...]) 和 BLAST+ 工具
AttributeError: 'PairwiseAligner' object has no attribute 'align'	Biopython < 1.78	升级：pip install --upgrade biopython（PairwiseAligner 要求 ≥ 1.72；从 1.78 起稳定）
Bio.pairwise2 弃用警告	使用旧版 pairwise2 API	替换为 from Bio.Align import PairwiseAligner（在 Biopython 1.84+ 中已移除）
ValueError: Sequence contains letters not in the alphabet	序列在调用 translate() 前包含非标准字符（例如 N、模糊码）	去除或替换模糊碱基；使用 seq.translate(table=1)，它可以处理模糊密码子
BLAST 结果为空 / NCBIXML.read() 引发 StopIteration	XML 为空或格式错误；网络超时；查询过短	检查查询序列长度（建议 >10 aa）；将 .read() 改为 .parse()，并检查 blast_record.alignments 的长度
Phylo.draw() 卡住或生成空白图	缺少 matplotlib 后端	在导入 Phylo 前调用 import matplotlib; matplotlib.use("Agg")，以用于无头环境
TreeConstruction 距离矩阵维度不匹配	比对中包含重复 ID	构建比对前对 SeqRecord ID 去重：{r.id: r for r in records}.values()

相关技能

• biopython-molecular-biology — 限制性消化、PCR 引物设计、蛋白质结构分析（Bio.PDB）、分子量和 Tm 计算；可与此技能互补
• pysam-genomic-files — SAM/BAM/CRAM 比对文件访问、pileup 和区域查询；用于处理读段比对
• scikit-bio — 生态多样性统计（UniFrac、Bray-Curtis）、排序分析（PCoA）和微生物组分析；在多样性分析方面比 Biopython 的系统发育功能更专业
• gget-genomic-databases — 无需设置 Entrez 工作流，即可一行查询基因、运行 BLAST 和查询数据库
• etetoolkit — 高级系统发育树可视化、使用 NCBI 分类法进行注释，以及生成出版级树图

参考资料

• Biopython 教程和 Cookbook — 全面的官方教程（第 2、5、6、7、9、11 章）
• Biopython API 参考 — 完整的模块和类文档
• GitHub：biopython/biopython — 源代码、变更日志和问题跟踪器
• NCBI E-utilities 文档 — esearch、efetch、elink、esummary 的 Entrez API 参考
• PairwiseAligner 迁移指南 — 如何替换已弃用的 Bio.pairwise2 模块