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() を使って大きな multi-FASTA または FASTQ ファイルをインデックス化し、すべての配列を RAM に読み込まずにランダムアクセスで取得する
• 1 回の呼び出しで配列形式を相互変換する（FASTA ↔ GenBank ↔ FASTQ ↔ PHYLIP）
• Newick または Nexus 系統樹をプログラムからたどり、ルート設定、枝刈り、注釈付けを行う
• SAM/BAM/CRAM アラインメントファイルやマッピング済みリードを扱う場合は、代わりに pysam を使用する
• 生態学的多様性指標（UniFrac、ベータ多様性）や、PCoA などの序列化手法が必要な場合は、代わりに scikit-bio を使用する
• Entrez のパイプラインを書かずに、遺伝子を手早く検索したり、NCBI に対するクエリを 1 行で実行したりする場合は、代わりに gget を使用する

前提条件

• Python パッケージ：biopython、numpy、matplotlib
• オプションのツール：MUSCLE または ClustalW をローカルにインストール（Bio.Align.Applications ラッパー用）
• NCBI へのアクセス：E-utilities を呼び出す前に Entrez.email を設定します。10 req/s（デフォルトは 3 req/s）で利用するには、https://www.ncbi.nlm.nih.gov/account/ で無料の API キーを取得してください
• ローカル BLAST：オフライン検索には BLAST+ を別途インストールします（conda install -c bioconda blast）

インストール前に確認してください：現在の環境（例：pixi / conda 環境内）でツールをすでに利用できる場合があります。まず command -v python を実行し、パスが返された場合は以下のインストールコマンドを省略してください。pixi プロジェクト内で実行する場合は、python 単体ではなく pixi run 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 req/s のレート制限（API キーを使用する場合は 10 req/s）を守ってください。

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}")

Module 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")

Module 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 オブジェクトのリストを格納します。フィーチャーでは、Python のスライスと同じ規約である、0 始まりの半開区間 FeatureLocation(start, end, strand) を使用します。鎖の向きを正しく反映した部分配列を取得するには、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 が設定されていることがあります。terminal=True/False を指定して tree.find_clades() を使用すると、対象を絞り込めます。ルートは tree.root（これも Clade です）です。Phylo.read() は Tree ラッパーを返し、tree.root からルートの Clade を取得できます。

一般的なワークフロー

ワークフロー 1：遺伝子ファミリーの系統解析 — NCBI からのダウンロード → MUSCLE アラインメント → 系統樹の構築

目標：Vertebrata の 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 ファイルを読み込み、平均品質と長さでリードをフィルタリングし、後続のアラインメントに使用する 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	負の浮動小数点数	ギャップを 1 残基伸長するごとのペナルティ
method	DistanceTreeConstructor	"nj"	"nj"、"upgma"	NJ（無根）と UPGMA（根付き、分子時計を仮定）の切り替え
format	Phylo.read/write()	必須	"newick"、"nexus"、"phyloxml"	系統樹ファイルの形式
rettype	Entrez.efetch()	必須	"fasta"、"gb"、"xml"	返されるレコードの形式。対応する SeqIO 形式文字列と組み合わせます

ベストプラクティス

1. 必ず Entrez.email を設定し、レート制限を守ってください。API キーなしで 3 req/s を超えると、NCBI によって IP アドレスがブロックされます。ループ内ではリクエストの間に time.sleep(0.4) を追加するか、Entrez.api_key を使って 10 req/s に対応してください。

   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 ファイルには list(SeqIO.parse(...)) ではなく、SeqIO.index() を使ってください。すべてのレコードをリストに読み込むと 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. 旧式の pairwise2 ではなく、PairwiseAligner を使ってください。Bio.pairwise2 は Biopython 1.80 以降で非推奨となっており、今後削除される予定です。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 経由で -outfmt 5（XML）を指定して blastp／blastn を呼び出し、NCBIXML で解析してください。

6. 距離を測定する前に系統樹をルート化してください。Phylo.distance() は、2 つのノード間の経路に沿った枝長の合計を測定します。ルート化されていない系統樹でも経路は一意ですが、可視化やクレードの判定には、中点ルート化または外群ルート化が生物学的に妥当です。

7. SeqRecord のコレクションを作成する前に、アラインメントのギャップを取り除いてください。BLAST やアラインメントの結果には、ギャップ文字（-）が含まれることがあります。.translate() などの Seq 操作は、ギャップを含む配列でエラーになります。seq.replace("-", "") で取り除くか、ungap() を使ってください。

よく使うレシピ

レシピ：履歴サーバーを使った 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 req/s）を超過しています	呼び出しの間に time.sleep(0.4) を追加するか、https://www.ncbi.nlm.nih.gov/account/ で無料の API キーを登録して 10 req/s にします
RuntimeError: Too many requests were made without pausing	NCBI が短時間に連続したリクエストを検出しています	Entrez.email を設定し（必須）、ループの頻度を下げ、ID をカンマ区切りの文字列にまとめて 1 回の 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 パイプラインを構築せずに、1 行で遺伝子検索、BLAST、データベースクエリを実行します
• etetoolkit — 高度な系統樹の可視化、NCBI 分類に基づく注釈、出版品質の系統樹の描画

参考資料

• Biopython チュートリアルとクックブック — 公式の包括的なチュートリアル（第 2、5、6、7、9、11 章）
• Biopython API リファレンス — モジュールとクラスの完全なドキュメント
• GitHub：biopython/biopython — ソース、変更履歴、課題トラッカー
• NCBI E-utilities ドキュメント — esearch、efetch、elink、esummary の Entrez API リファレンス
• PairwiseAligner 移行ガイド — 非推奨の Bio.pairwise2 モジュールの置き換え