biopython-molecular-biology について
Biopythonは、計算分子生物学向けのオープンソースPythonライブラリです。生物配列、データベース、構造、系統解析を扱う研究者や開発者を対象とし、ファイル解析、NCBIアクセス、BLAST、アラインメントなどのAPIを備えています。
主な機能
- 01生物学的ファイル形式の解析と変換
- 02NCBIデータベースへのプログラムアクセス
- 03BLAST検索の実行と解析
- 04ペアワイズ配列アラインメント
- 05タンパク質構造の解析
- 06系統樹の構築と可視化
活用シーン
- 01遺伝子配列の取得と分析
- 02配列比較と系統樹の作成
- 03大規模なFASTQまたはFASTAデータの選別と集計
- 04タンパク質構造の確認とRMSD計算
インストール
コマンドをコピーしてターミナルで実行します。
npx skills add https:/ / github.com/ jaechang-hits/ sciagent-skills --skill biopython-molecular-biologyインストール後、Claude Code で/biopython-molecular-biologyを入力して読み込まれたか確認します。
詳細
Biopython:計算分子生物学ツールキット
概要
Biopython は、計算分子生物学向けの標準的なオープンソース Python ライブラリです。配列の処理、生物学的ファイルの解析、NCBI データベースへのアクセス、BLAST 検索、タンパク質構造解析、系統学のためのモジュール式 API を提供します。Python 3 をサポートし、NumPy が必要です。
使用する場面
- 生物学的ファイル形式(FASTA、GenBank、FASTQ、PDB、mmCIF、PHYLIP)を解析・変換する
- NCBI データベース(GenBank、PubMed、Protein)から配列や文献をプログラムで取得する
- BLAST 検索(NCBI のリモートサーバーまたはローカルの BLAST+)を実行・解析する
- 独自のスコアリングを使って、ペアワイズまたはマルチプル配列アラインメントを実行する
- 3D タンパク質構造を解析する — 距離、角度、DSSP、重ね合わせ
- 配列アラインメントから系統樹を構築・可視化する
- 配列の統計値(GC 含量、分子量、融解温度)を計算する
- 独自のフィルタリングロジックを使って、数千の配列をバッチ処理する
- SAM/BAM/CRAM アラインメントファイルの読み込みやマッピング済みリードの処理には、代わりに
pysamを使用します。高度な生態学的多様性指標には、代わりにscikit-bioを使用します
前提条件
- Python パッケージ:
biopython、numpy(樹形図の可視化にはmatplotlib) - データ要件:配列ファイル(FASTA、GenBank、FASTQ)または NCBI アクセス用のアクセッション ID
- 環境:Python 3.8 以降。NCBI Entrez の利用にはメールアドレスの登録が必要です
pip install biopython numpy matplotlib
クイックスタート
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction
# Parse a FASTA file and compute basic statistics
records = list(SeqIO.parse("sequences.fasta", "fasta"))
print(f"Sequences loaded: {len(records)}")
seq = records[0].seq
print(f"ID: {records[0].id}")
print(f"Length: {len(seq)} bp")
print(f"GC content: {gc_fraction(seq)*100:.1f}%")
print(f"Reverse complement: {seq.reverse_complement()[:30]}...")
print(f"Protein translation: {seq.translate()[:10]}...")
コア API
モジュール 1:配列オブジェクト(Bio.Seq)
DNA、RNA、タンパク質の配列を作成・操作します。
from Bio.Seq import Seq
# Create sequence and perform standard operations
dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
print(f"Length: {len(dna)} bp")
print(f"Complement: {dna.complement()}")
print(f"Reverse complement: {dna.reverse_complement()}")
print(f"Transcription: {dna.transcribe()}")
print(f"Translation: {dna.translate()}")
print(f"Translation (to stop): {dna.translate(to_stop=True)}")
# Length: 39 bp
# Translation: MAIVMGR*KGAR*
# Translation (to stop): MAIVMGR
from Bio.Seq import Seq
# Alternative genetic codes (e.g., mitochondrial)
mito_dna = Seq("ATGGCCATTGTAATGGGCCGCTGA")
std_protein = mito_dna.translate(table=1) # Standard
mito_protein = mito_dna.translate(table=2) # Vertebrate mitochondrial
print(f"Standard: {std_protein}")
print(f"Mitochondrial: {mito_protein}")
# Find all start codons
coding_dna = Seq("ATGAAACCCATGGGGTTTAAATAG")
positions = [i for i in range(len(coding_dna) - 2) if coding_dna[i:i+3] == "ATG"]
print(f"ATG positions: {positions}")
# ATG positions: [0, 9]
モジュール 2:配列 I/O(Bio.SeqIO)
生物学的ファイル形式を読み込み、書き込み、変換します。
from Bio import SeqIO
# Parse FASTA file — returns SeqRecord iterator
records = list(SeqIO.parse("sequences.fasta", "fasta"))
print(f"Loaded {len(records)} sequences")
for rec in records[:3]:
print(f" {rec.id}: {len(rec.seq)} bp — {rec.description}")
# Parse GenBank — rich annotation access
for rec in SeqIO.parse("genome.gb", "genbank"):
print(f"{rec.id}: {len(rec.features)} features, {len(rec.seq)} bp")
for feat in rec.features[:5]:
print(f" {feat.type}: {feat.location}")
# Convert between formats
count = SeqIO.convert("input.gb", "genbank", "output.fasta", "fasta")
print(f"Converted {count} records: GenBank → FASTA")
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq
# Write sequences to file
records = [
SeqRecord(Seq("ATCGATCG"), id="seq1", description="test sequence 1"),
SeqRecord(Seq("GCTAGCTA"), id="seq2", description="test sequence 2"),
]
count = SeqIO.write(records, "output.fasta", "fasta")
print(f"Wrote {count} records to output.fasta")
# Filter sequences by length (streaming — memory efficient)
long_seqs = (rec for rec in SeqIO.parse("large_file.fasta", "fasta") if len(rec.seq) >= 200)
count = SeqIO.write(long_seqs, "filtered.fasta", "fasta")
print(f"Kept {count} sequences >= 200 bp")
# Index large FASTA for random access
idx = SeqIO.index("large_file.fasta", "fasta")
print(f"Indexed {len(idx)} sequences")
rec = idx["target_sequence_id"]
print(f"Retrieved: {rec.id}, {len(rec.seq)} bp")
モジュール 3:NCBI データベースへのアクセス(Bio.Entrez)
NCBI データベースをプログラムで検索し、データをダウンロードします。
from Bio import Entrez, SeqIO
Entrez.email = "[email protected]"
# Entrez.api_key = "your_key" # Optional: 10 req/s instead of 3 req/s
# Search PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR Cas9 2024", retmax=5)
results = Entrez.read(handle)
handle.close()
print(f"Found {results['Count']} articles, retrieved {len(results['IdList'])} IDs")
print(f"IDs: {results['IdList']}")
# Fetch GenBank record by accession
handle = Entrez.efetch(db="nucleotide", id="EU490707", rettype="gb", retmode="text")
record = SeqIO.read(handle, "genbank")
handle.close()
print(f"{record.id}: {record.description}")
print(f"Length: {len(record.seq)} bp, Features: {len(record.features)}")
from Bio import Entrez
import time
Entrez.email = "[email protected]"
# Batch download with rate limiting
handle = Entrez.esearch(db="protein", term="insulin[Protein Name] AND human[Organism]", retmax=20)
results = Entrez.read(handle)
handle.close()
# Fetch summaries in batch
ids = results["IdList"][:10]
handle = Entrez.esummary(db="protein", id=",".join(ids))
summaries = Entrez.read(handle)
handle.close()
for doc in summaries:
print(f" {doc['AccessionVersion']}: {doc['Title'][:60]}...")
print(f"\nFetched {len(summaries)} protein summaries")
モジュール 4:BLAST 操作(Bio.Blast)
NCBI またはローカルのデータベースに対して BLAST 検索を実行・解析します。
from Bio.Blast import NCBIWWW, NCBIXML
# Remote BLAST search (nucleotide)
query_seq = "ATCGATCGATCGATCGATCGATCGATCGATCG"
result_handle = NCBIWWW.qblast("blastn", "nt", query_seq, hitlist_size=5)
blast_record = NCBIXML.read(result_handle)
result_handle.close()
print(f"Query: {blast_record.query[:50]}")
print(f"Database: {blast_record.database}")
print(f"Hits: {len(blast_record.alignments)}")
for aln in blast_record.alignments[:3]:
hsp = aln.hsps[0]
print(f"\n {aln.title[:60]}...")
print(f" E-value: {hsp.expect:.2e}, Identity: {hsp.identities}/{hsp.align_length}")
print(f" Score: {hsp.score}, Bits: {hsp.bits:.1f}")
from Bio.Blast.Applications import NcbiblastpCommandline
from Bio.Blast import NCBIXML
# Local BLAST (requires BLAST+ installed)
blastp_cline = NcbiblastpCommandline(
query="query.fasta",
db="swissprot",
evalue=1e-5,
outfmt=5, # XML output
out="blast_results.xml",
num_threads=4,
)
print(f"Command: {blastp_cline}")
# stdout, stderr = blastp_cline() # Execute
# Parse local BLAST XML results
with open("blast_results.xml") as f:
for record in NCBIXML.parse(f):
print(f"Query: {record.query}")
for aln in record.alignments[:3]:
print(f" Hit: {aln.hit_def[:50]}, E={aln.hsps[0].expect:.2e}")
モジュール 5:ペアワイズアラインメント(Bio.Align)
スコアリングをカスタマイズして、グローバルまたはローカルのペアワイズ配列アラインメントを実行します。
from Bio import Align
# Global alignment
aligner = Align.PairwiseAligner()
aligner.mode = "global"
aligner.match_score = 2
aligner.mismatch_score = -1
aligner.open_gap_score = -5
aligner.extend_gap_score = -0.5
alignments = aligner.align("ACCGGTAACG", "ACGGTAAC")
print(f"Score: {alignments.score}")
print(f"Number of alignments: {len(alignments)}")
print(f"Best alignment:\n{alignments[0]}")
from Bio import Align
from Bio.Align import substitution_matrices
# Protein alignment with BLOSUM62
aligner = Align.PairwiseAligner()
aligner.mode = "local"
aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
aligner.open_gap_score = -10
aligner.extend_gap_score = -0.5
seq1 = "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH"
seq2 = "MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLST"
alignments = aligner.align(seq1, seq2)
print(f"Score: {alignments.score}")
print(f"Best alignment:\n{alignments[0]}")
モジュール 6:タンパク質構造解析(Bio.PDB)
PDB/mmCIF ファイルを解析し、3D タンパク質構造を分析します。
from Bio.PDB import PDBParser, PPBuilder
# Parse PDB structure
parser = PDBParser(QUIET=True)
structure = parser.get_structure("1CRN", "1crn.pdb")
# Navigate SMCRA hierarchy: Structure > Model > Chain > Residue > Atom
model = structure[0]
for chain in model:
residues = list(chain.get_residues())
print(f"Chain {chain.id}: {len(residues)} residues")
# Extract sequence from structure
ppb = PPBuilder()
for pp in ppb.build_peptides(structure):
print(f"Peptide: {pp.get_sequence()[:50]}... ({len(pp.get_sequence())} aa)")
# Calculate CA-CA distance
chain_a = model["A"]
ca1 = chain_a[10]["CA"]
ca2 = chain_a[20]["CA"]
distance = ca1 - ca2
print(f"CA distance (res 10-20): {distance:.2f} Angstrom")
from Bio.PDB import PDBParser, Superimposer
import numpy as np
# Structure superimposition (RMSD calculation)
parser = PDBParser(QUIET=True)
struct1 = parser.get_structure("s1", "structure1.pdb")
struct2 = parser.get_structure("s2", "structure2.pdb")
# Get CA atoms for alignment
atoms1 = [res["CA"] for res in struct1[0]["A"].get_residues() if "CA" in res]
atoms2 = [res["CA"] for res in struct2[0]["A"].get_residues() if "CA" in res]
# Superimpose (requires same number of atoms)
n = min(len(atoms1), len(atoms2))
sup = Superimposer()
sup.set_atoms(atoms1[:n], atoms2[:n])
sup.apply(struct2.get_atoms())
print(f"RMSD: {sup.rms:.3f} Angstrom over {n} CA atoms")
モジュール 7:系統学(Bio.Phylo)
系統樹を構築、操作、可視化します。
from Bio import Phylo, AlignIO
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
import io
# Build tree from multiple sequence alignment
alignment = AlignIO.read("aligned_sequences.fasta", "fasta")
print(f"Alignment: {len(alignment)} sequences, {alignment.get_alignment_length()} positions")
# Calculate distance matrix and build NJ tree
calculator = DistanceCalculator("identity")
dm = calculator.get_distance(alignment)
print(f"Distance matrix:\n{dm}")
constructor = DistanceTreeConstructor()
nj_tree = constructor.nj(dm)
upgma_tree = constructor.upgma(dm)
# Visualize
Phylo.draw_ascii(nj_tree)
# Save tree
Phylo.write(nj_tree, "tree.nwk", "newick")
print("Saved tree.nwk")
モジュール 8:配列ユーティリティ(Bio.SeqUtils)
配列の統計値と物理化学的特性を計算します。
from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction, molecular_weight, MeltingTemp
dna = Seq("ATCGATCGATCGATCGATCG")
print(f"GC content: {gc_fraction(dna):.2%}")
print(f"Molecular weight: {molecular_weight(dna, seq_type='DNA'):.2f} Da")
print(f"Melting temp (basic): {MeltingTemp.Tm_Wallace(dna):.1f} C")
print(f"Melting temp (NN): {MeltingTemp.Tm_NN(dna):.1f} C")
# Protein analysis
from Bio.SeqUtils.ProtParam import ProteinAnalysis
protein = ProteinAnalysis("MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTK")
print(f"\nProtein MW: {protein.molecular_weight():.2f} Da")
print(f"Isoelectric point: {protein.isoelectric_point():.2f}")
print(f"Aromaticity: {protein.aromaticity():.4f}")
print(f"Instability index: {protein.instability_index():.2f}")
print(f"GRAVY: {protein.gravy():.4f}")
aa_pct = protein.get_amino_acids_percent()
print(f"Top 3 amino acids: {sorted(aa_pct.items(), key=lambda x: -x[1])[:3]}")
一般的なワークフロー
ワークフロー 1:遺伝子配列の取得と解析のパイプライン
目的:NCBI から遺伝子を取得し、その特性を解析してタンパク質に翻訳し、統計値を計算します。
from Bio import Entrez, SeqIO
from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction
from Bio.SeqUtils.ProtParam import ProteinAnalysis
Entrez.email = "[email protected]"
# Step 1: Fetch gene from GenBank
handle = Entrez.efetch(db="nucleotide", id="NM_007294.4", rettype="gb", retmode="text")
record = SeqIO.read(handle, "genbank")
handle.close()
print(f"Gene: {record.description}")
print(f"Length: {len(record.seq)} bp")
# Step 2: Extract CDS
cds_features = [f for f in record.features if f.type == "CDS"]
if cds_features:
cds = cds_features[0]
cds_seq = cds.location.extract(record).seq
print(f"CDS: {len(cds_seq)} bp, GC: {gc_fraction(cds_seq):.2%}")
# Step 3: Translate
protein_seq = cds_seq.translate(to_stop=True)
print(f"Protein: {len(protein_seq)} aa")
print(f"First 30 aa: {protein_seq[:30]}...")
# Step 4: Protein properties
analysis = ProteinAnalysis(str(protein_seq))
print(f"MW: {analysis.molecular_weight():.0f} Da")
print(f"pI: {analysis.isoelectric_point():.2f}")
print(f"Instability: {analysis.instability_index():.1f}")
print(f"GRAVY: {analysis.gravy():.3f}")
ワークフロー 2:BLAST と系統解析による配列の比較分析
目的:タンパク質配列を BLAST 検索し、相同配列を取得してアラインメントを行い、系統樹を構築します。
from Bio.Blast import NCBIWWW, NCBIXML
from Bio import Entrez, SeqIO, AlignIO, Phylo
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Align.Applications import MuscleCommandline
import time
Entrez.email = "[email protected]"
# Step 1: BLAST search
query = "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH"
result_handle = NCBIWWW.qblast("blastp", "swissprot", query, hitlist_size=10)
blast_record = NCBIXML.read(result_handle)
print(f"BLAST hits: {len(blast_record.alignments)}")
# Step 2: Collect homolog accessions
accessions = []
for aln in blast_record.alignments[:8]:
acc = aln.accession
accessions.append(acc)
print(f" {acc}: E={aln.hsps[0].expect:.2e}, {aln.hit_def[:50]}...")
# Step 3: Fetch sequences and save for alignment
handle = Entrez.efetch(db="protein", id=",".join(accessions), rettype="fasta", retmode="text")
records = list(SeqIO.parse(handle, "fasta"))
handle.close()
SeqIO.write(records, "homologs.fasta", "fasta")
print(f"Saved {len(records)} homolog sequences")
# Step 4: Align (requires MUSCLE installed)
# muscle_cline = MuscleCommandline(input="homologs.fasta", out="aligned.fasta")
# muscle_cline()
# Step 5: Build phylogenetic tree from alignment
alignment = AlignIO.read("aligned.fasta", "fasta")
calculator = DistanceCalculator("blosum62")
dm = calculator.get_distance(alignment)
tree = DistanceTreeConstructor().nj(dm)
Phylo.draw_ascii(tree)
Phylo.write(tree, "homologs.nwk", "newick")
print("Saved phylogenetic tree to homologs.nwk")
ワークフロー 3:配列のバッチ処理と品質フィルタリング
目的:大規模な FASTQ/FASTA データセットを処理し、品質や長さでフィルタリングして統計値を計算し、エクスポートします。
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction
import pandas as pd
# Step 1: Stream through large file and collect stats
stats = []
passed = []
for rec in SeqIO.parse("reads.fastq", "fastq"):
seq_len = len(rec.seq)
gc = gc_fraction(rec.seq)
avg_qual = sum(rec.letter_annotations["phred_quality"]) / seq_len
stats.append({"id": rec.id, "length": seq_len, "gc": gc, "avg_qual": avg_qual})
# Filter: length >= 100 and avg quality >= 20
if seq_len >= 100 and avg_qual >= 20:
passed.append(rec)
# Step 2: Summary statistics
df = pd.DataFrame(stats)
print(f"Total reads: {len(df)}")
print(f"Passed QC: {len(passed)} ({len(passed)/len(df)*100:.1f}%)")
print(f"\nLength: mean={df['length'].mean():.0f}, median={df['length'].median():.0f}")
print(f"GC: mean={df['gc'].mean():.2%}, std={df['gc'].std():.2%}")
print(f"Quality: mean={df['avg_qual'].mean():.1f}, min={df['avg_qual'].min():.1f}")
# Step 3: Export filtered reads
count = SeqIO.write(passed, "filtered_reads.fastq", "fastq")
print(f"\nExported {count} filtered reads to filtered_reads.fastq")
# Step 4: Save statistics
df.to_csv("read_statistics.csv", index=False)
print(f"Saved statistics to read_statistics.csv")
主なパラメーター
| パラメーター | モジュール | デフォルト | 範囲/選択肢 | 効果 |
|---|---|---|---|---|
Seq.translate(table=) | Bio.Seq | 1(標準) | 1~33 | 翻訳に使用するNCBI遺伝暗号表 |
Seq.translate(to_stop=) | Bio.Seq | False | True、False | 最初の終止コドンで停止するか、配列全体を翻訳するか |
SeqIO.parse(format=) | Bio.SeqIO | — | "fasta"、"genbank"、"fastq"、"phylip" | パースするファイル形式 |
Entrez.email | Bio.Entrez | (必須) | 有効なメールアドレス | NCBIは追跡のためにメールアドレスを必要とします。Entrezを呼び出す前に設定してください |
Entrez.api_key | Bio.Entrez | None | NCBI APIキー文字列 | リクエストレート制限が毎秒3件から10件に増加します |
PairwiseAligner.mode | Bio.Align | "global" | "global"、"local" | Needleman-WunschアルゴリズムとSmith-Watermanアルゴリズムのどちらを使用するか |
PairwiseAligner.open_gap_score | Bio.Align | -1 | -20~0 | ギャップ開始時のペナルティ。値がより負であるほどギャップが少なくなります |
PairwiseAligner.substitution_matrix | Bio.Align | None | "BLOSUM62"、"BLOSUM45"、"PAM250" | タンパク質アラインメント用のスコアリング行列 |
PDBParser(QUIET=) | Bio.PDB | False | True、False | 標準外のPDBファイルに関するパーサーの警告を抑制します |
DistanceCalculator(model=) | Bio.Phylo | "identity" | "identity"、"blosum62" | 系統樹構築に使用する距離モデル |
よく使うレシピ
レシピ:制限酵素解析
使用する場面:クローニング設計のために、DNA配列内の制限酵素認識部位を検索します。
from Bio.Seq import Seq
from Bio.Restriction import EcoRI, BamHI, HindIII, RestrictionBatch, Analysis
seq = Seq("GAATTCAAAGGATCCTTTTAAGCTTGGGAATTC")
# Single enzyme
print(f"EcoRI cuts at: {EcoRI.search(seq)}")
print(f"BamHI cuts at: {BamHI.search(seq)}")
# Batch analysis
batch = RestrictionBatch([EcoRI, BamHI, HindIII])
analysis = Analysis(batch, seq)
result = analysis.full()
for enzyme, sites in result.items():
if sites:
print(f"{enzyme}: cuts at positions {sites}")
レシピ:モチーフ探索と位置重み行列
使用する場面:転写因子結合部位やコンセンサスパターンを解析します。
from Bio import motifs
from Bio.Seq import Seq
# Create motif from observed binding sites
instances = [
Seq("TACGAT"),
Seq("TAGCAT"),
Seq("TACGGT"),
Seq("TAGCAT"),
Seq("TACGAT"),
]
m = motifs.create(instances)
print(f"Consensus: {m.consensus}")
print(f"Degenerate: {m.degenerate_consensus}")
print(f"\nPosition Weight Matrix:")
print(m.counts)
# Score a new sequence against the motif
pwm = m.counts.normalize(pseudocounts=0.5)
pssm = pwm.log_odds()
test_seq = Seq("AATACGATCCC")
for pos, score in pssm.search(test_seq, threshold=0.0):
print(f" Position {pos}: score={score:.2f}")
レシピ:GenomeDiagram — ゲノム上の特徴を可視化
使用する場面:GenBankファイルから円形または線形のゲノムマップを作成します。
from Bio import SeqIO
from Bio.Graphics import GenomeDiagram
from reportlab.lib import colors
from reportlab.lib.units import cm
# Load annotated genome
record = SeqIO.read("plasmid.gb", "genbank")
# Create diagram
diagram = GenomeDiagram.Diagram(record.name)
track = diagram.new_track(1, name="Annotated Features", greytrack=True)
feature_set = track.new_set()
# Color features by type
color_map = {"CDS": colors.blue, "gene": colors.green, "promoter": colors.red}
for feature in record.features:
if feature.type in color_map:
feature_set.add_feature(feature, color=color_map[feature.type], label=True, label_size=8)
# Draw circular map
diagram.draw(format="circular", circular=True, pagesize=(20*cm, 20*cm), start=0, end=len(record))
diagram.write("genome_map.pdf", "PDF")
print(f"Saved genome_map.pdf ({len(record.features)} features)")
レシピ:PubMed文献の一括検索
使用する場面:プログラムから文献メタデータを検索してダウンロードします。
from Bio import Entrez
import time
Entrez.email = "[email protected]"
# Search PubMed with complex query
query = "(CRISPR[Title]) AND (2024[Date - Publication]) AND (review[Publication Type])"
handle = Entrez.esearch(db="pubmed", term=query, retmax=100)
results = Entrez.read(handle)
handle.close()
print(f"Found {results['Count']} articles")
# Fetch abstracts in batches
ids = results["IdList"]
batch_size = 20
articles = []
for i in range(0, len(ids), batch_size):
batch = ids[i:i+batch_size]
handle = Entrez.efetch(db="pubmed", id=",".join(batch), rettype="xml")
records = Entrez.read(handle)
handle.close()
for article in records["PubmedArticle"]:
info = article["MedlineCitation"]["Article"]
title = info.get("ArticleTitle", "N/A")
abstract = info.get("Abstract", {}).get("AbstractText", ["N/A"])[0]
articles.append({"title": title, "abstract": str(abstract)[:200]})
time.sleep(0.4) # Rate limiting
for a in articles[:5]:
print(f" {a['title'][:70]}...")
print(f"\nRetrieved {len(articles)} articles")
トラブルシューティング
| 問題 | 原因 | 解決策 |
|---|---|---|
EntrezからのHTTPError 400 | アクセッション/IDが無効、またはクエリの形式が不正 | アクセッションを検証します。まずNCBIのウェブインターフェースでクエリ構文を確認してください |
HTTPError 429(リクエスト過多) | NCBIのレート制限(毎秒3件)を超過 | Entrez.api_keyを設定して毎秒10件にします。呼び出しの間にtime.sleep(0.4)を追加してください |
SeqIO.read()でValueError: No records found | ファイルが空、または形式文字列が誤っている | SeqIO.parse()を使ってファイルにレコードがあるか確認します。形式が実際の内容と一致しているか確認してください |
Bio.PDB.PDBExceptions.PDBConstructionWarning | 標準外の原子、または占有率の問題 | PDBParser(QUIET=True)を使うか、pdb-toolsでPDBを修正します。代替コンフォメーションがないか確認してください |
| BLAST検索がタイムアウトする | クエリが大きい、またはNCBIサーバーが混雑している | 大規模な検索にはローカルのBLAST+を使います。NCBIWWW.qblast(hitlist_size=N)を設定して結果数を制限してください |
Alignment has sequences of different lengths | アラインメントされていない配列をAlignIOに渡した | 配列をアラインメントとして読み込む前に、MUSCLE/Clustalでアラインメントしてください |
SeqIO.index()がValueErrorを発生させる | FASTAファイルに重複IDがある | SeqIO.to_dict()を使ってIDを重複排除するか、awkで前処理してください |
ImportError: No module named Bio | 有効な環境にBiopythonがインストールされていない | pip install biopythonを実行します。python -c "import Bio; print(Bio.__version__)"で確認してください |
translate()の結果に予期しない*が含まれる | 配列の途中に終止コドンがある | 読み枠を確認します。正しい遺伝暗号を使ってSeq.translate(table=N)を指定してください |
| GenomeDiagramの出力が空白 | フィルター条件に一致する特徴がない | GenBankファイル内のfeature.typeの値を確認します。デバッグのために型を出力してください |
参考資料
- Biopython Tutorial and Cookbook — 公式の包括的なチュートリアル
- Biopython API Documentation — API の完全なリファレンス
- Cock, P.J.A. et al. (2009) Bioinformatics 25:1422-1423 — Biopython に関する原著論文
- NCBI Entrez Programming Utilities — NCBI へのプログラムによるアクセスに関する E-utilities のドキュメント
- Biopython GitHub Repository — ソースコード、イシュー、リリースノート