mcpField
B
biopython-molecular-biology

biopython-molecular-biology

github.com/jaechang-hitsデータ・データベース38482 回インストール6 日前に更新
リポジトリ

biopython-molecular-biology について

Biopythonは、計算分子生物学向けのオープンソースPythonライブラリです。生物配列、データベース、構造、系統解析を扱う研究者や開発者を対象とし、ファイル解析、NCBIアクセス、BLAST、アラインメントなどのAPIを備えています。

主な機能

  1. 01生物学的ファイル形式の解析と変換
  2. 02NCBIデータベースへのプログラムアクセス
  3. 03BLAST検索の実行と解析
  4. 04ペアワイズ配列アラインメント
  5. 05タンパク質構造の解析
  6. 06系統樹の構築と可視化

活用シーン

  1. 01遺伝子配列の取得と分析
  2. 02配列比較と系統樹の作成
  3. 03大規模なFASTQまたはFASTAデータの選別と集計
  4. 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.Seq1(標準)1~33翻訳に使用するNCBI遺伝暗号表
Seq.translate(to_stop=)Bio.SeqFalseTrue、False最初の終止コドンで停止するか、配列全体を翻訳するか
SeqIO.parse(format=)Bio.SeqIO—"fasta"、"genbank"、"fastq"、"phylip"パースするファイル形式
Entrez.emailBio.Entrez(必須)有効なメールアドレスNCBIは追跡のためにメールアドレスを必要とします。Entrezを呼び出す前に設定してください
Entrez.api_keyBio.EntrezNoneNCBI APIキー文字列リクエストレート制限が毎秒3件から10件に増加します
PairwiseAligner.modeBio.Align"global""global"、"local"Needleman-WunschアルゴリズムとSmith-Watermanアルゴリズムのどちらを使用するか
PairwiseAligner.open_gap_scoreBio.Align-1-20~0ギャップ開始時のペナルティ。値がより負であるほどギャップが少なくなります
PairwiseAligner.substitution_matrixBio.AlignNone"BLOSUM62"、"BLOSUM45"、"PAM250"タンパク質アラインメント用のスコアリング行列
PDBParser(QUIET=)Bio.PDBFalseTrue、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の値を確認します。デバッグのために型を出力してください

参考資料

コメント