忍者ブログ
バイオインフォマティックス技術者試験、情報処理試験など、IT系の試験を基礎から勉強します。また、Javaなどプログラミングを勉強します。

【バイオインフォ実習】第5回:NCBI DBから複数配列を自動検索・一括取得!マルチFASTAの活用とSeqIO.parseによるループ統計解析

前回(第4回)は、あらかじめ判明している1つのアクセッション番号(TP53)を指定して、NCBI APIから単一データを取得する方法を解説しました。

しかし、実際のデータ解析の現場では、「特定の検索条件にヒットする複数の配列を一括取得し、まとめて統計解析や比較解析を行う」というフローが頻出します。

今回は、NCBIデータベースへのキーワード検索(esearch)から、複数データのマルチFASTA一括ダウンロード(efetch)、そしてBiopythonの SeqIO.parse() を使ったループ統計処理までを完全自動化するパイプラインを構築します。

1. マルチFASTA(Multi-FASTA)とは?

マルチFASTAとは、1つのファイル内に複数のDNA・RNA・タンパク質配列を順番に格納したフォーマットです。

>NC_013993.1 Homo sp. Altai mitochondrion, complete genome
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT...
>NC_012920.1 Homo sapiens mitochondrion, complete genome
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT...
>NC_011137.1 Homo sapiens neanderthalensis mitochondrion, complete genome
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT...

単一FASTAとマルチFASTAの比較

  • 単一FASTA: 1ファイル = 1配列SeqIO.read() で読み込む)
  • マルチFASTA: 1ファイル = 複数の配列SeqIO.parse() でループ処理する)

複数種の比較解析(マルチプルアライメントや系統樹作成)や、特定の遺伝子群を一括処理する際に必須となる形式です。

2. 検索から一括取得・ループ解析までの完全コード

以下のスクリプトは、NCBIの nuccore データベースから「ヒト属のミトコンドリア完全ゲノム(RefSeq)」を検索し、ヒットした上位3件をマルチFASTA形式で自動取得して解析するコードです。

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

# 1. 共通設定
Entrez.email = "your_email@example.com" # 自身のメールアドレスに変更してください

# 検索クエリ(ヒト属のRefSeqミトコンドリア完全ゲノム)
search_term = "Homo sapiens[Organism] AND mitochondrion[Title] AND srcdb_refseq[PROP]"
max_results = 3 # 取得件数上限

print(f"検索クエリ: {search_term}")
print("NCBI データベースを検索中...")

# 2. NCBI 検索(esearch): 条件に合う ID リストの取得
with Entrez.esearch(db="nuccore", term=search_term, retmax=max_results) as handle:
    search_results = Entrez.read(handle)

id_list = search_results["IdList"]
count = search_results["Count"]

print(f"検索ヒット件数: {count} 件")
print(f"取得対象 ID リスト (上位{len(id_list)}件): {id_list}\n")

if not id_list:
    print("該当するデータが見つかりませんでした。")
    exit()

# 3. NCBI 一括取得(efetch): 複数 ID を指定してマルチFASTAを受信
print("マルチFASTAデータを一括ダウンロード中...")

# id_list(GI番号等のリスト)をカンマ区切り文字列にして API に送る
with Entrez.efetch(db="nuccore", id=",".join(id_list), rettype="fasta", retmode="text") as handle:
    # 取得したマルチFASTAデータをローカルファイルに保存
    output_filename = "ncbi_multi_records.fasta"
    with open(output_filename, "w") as out_f:
        out_f.write(handle.read())

print(f"保存完了: {output_filename}\n")

# 4. SeqIO.parse() によるマルチFASTAのループ解析
print(f"=== {output_filename} の解析結果 ===")

# 複数配列が含まれるため SeqIO.parse を使用(イテレータ処理)
total_length = 0
record_count = 0

for record in SeqIO.parse(output_filename, "fasta"):
    record_count += 1
    seq_len = len(record.seq)
    total_length += seq_len
    
    # GC含有率の計算
    gc_val = gc_fraction(record.seq) * 100
    
    print(f"[{record_count}] ID: {record.id}")
    print(f" 概要: {record.description[:60]}...")
    print(f" 配列長: {seq_len:,} bp")
    print(f" GC率 : {gc_val:.2f}%")
    print("-" * 60)

# 全体の要約統計
avg_length = total_length / record_count if record_count > 0 else 0
print(f"解析完了: 合計 {record_count} 件 | 平均配列長: {avg_length:,.1f} bp")

3. 実行結果の確認

上記のプログラムを実行すると、アルタイ人(デニソワ人洞窟の個体)、現代人、ネアンデルタール人のミトコンドリア全ゲノムデータが自動で取得され、以下のようにパースされます。

検索クエリ: Homo sapiens[Organism] AND mitochondrion[Title] AND srcdb_refseq[PROP]
NCBI データベースを検索中...
検索ヒット件数: 3 件
取得対象 ID リスト (上位3件): ['292606408', '251831106', '196123578']

マルチFASTAデータを一括ダウンロード中...
保存完了: ncbi_multi_records.fasta

=== ncbi_multi_records.fasta の解析結果 ===
[1] ID: NC_013993.1
    概要: NC_013993.1 Homo sp. Altai mitochondrion, complete genome...
    配列長: 16,570 bp
    GC率 : 44.31%
------------------------------------------------------------
[2] ID: NC_012920.1
    概要: NC_012920.1 Homo sapiens mitochondrion, complete genome...
    配列長: 16,569 bp
    GC率 : 44.36%
------------------------------------------------------------
[3] ID: NC_011137.1
    概要: NC_011137.1 Homo sapiens neanderthalensis mitochondrion, complete...
    配列長: 16,565 bp
    GC率 : 44.39%
------------------------------------------------------------
解析完了: 合計 3 件 | 平均配列長: 16,568.0 bp

4. コードとデータ処理のポイント

一連のコードの中で、APIとBiopythonがどのように連動しているかをエンジニア目線で解説します。

  • esearch による内部ID(GI番号)の取得:
    Entrez.esearch を実行すると、search_results["IdList"]['292606408', '251831106', '196123578'] というGI番号(NCBI内部の識別ID)が動的に返されます。
  • カンマ区切りによる efetch の一括リクエスト:
    id=",".join(id_list) により、複数のGI番号を1つのリクエストとして送信しています。NCBIサーバーからまとめてマルチFASTA形式で返却されるため、通信回数を削減して効率的なダウンロードが可能です。
  • SeqIO.parse() によるメモリ効率の良いループ処理:
    SeqIO.read() は1ファイル1配列専用ですが、SeqIO.parse() は配列を1つずつ順番に読み込むイテレータ(Generator)として機能します。そのため、数百〜数千件の配列が含まれる大規模ファイルでもメモリを消費せずに高速処理できます。

5. ITエンジニア的まとめ

今回の実習のポイントは以下の通りです。

  • 自動検索と一括取得: esearch で得たIDリストをカンマ連結して efetch に送ることで、マルチFASTAの一括ダウンロードが可能。
  • イテレータ処理によるパイプライン化: 複数配列の解析には SeqIO.parse() を使用し、for ループで各レコードの塩基長やGC率を動的計算する。

次回は【バイオインフォ実習】第6回:GenBank形式ファイルからのアノテーション抽出とPandas DataFrame化に挑戦します!




PR

【バイオインフォ実習】第4回:Bio.EntrezでAPI連携!FASTAデータの自動取得パイプライン構築とRefSeq採番規則の完全理解

前回(第3回)は、FASTA形式の基本構造とBiopythonの SeqRecord オブジェクトによるパース方法を学びました。第4回の今回は、いよいよPythonコードからNCBI(米国国立生物工学情報センター)のAPIを呼び出し、本物の遺伝子データ(FASTA)を自動取得・解析するパイプラインを作成します。

Webブラウザでダウンロードする手間を省き、コード一発で公共データベースと連携する「ITエンジニアらしい実務的な手法」をマスターしていきましょう。

1. 今回扱うサンプルデータ:ヒトTP53遺伝子(NM_000546.6)

今回取得するのは、ヒトのTP53遺伝子(mRNA)のデータです。この NM_000546.6 という記号は、NCBIが管理するRefSeq(Reference Sequence)データベースにおける固有のルール(採番規則)に従って厳密に決められています。単なるランダムなIDではなく、コード自体に意味が埋め込まれています。

① 記号の解読(RefSeqの採番規則)

  • 1. プレフィックス NM_(mRNAの意味):
    先頭の2文字+アンダースコアは、データの種類(分子種)を表しています。
    • NM_:RefSeq curated mRNA(構造が手動・実験的に検証・整理された「mRNA」の配列データ)
    • NP_:RefSeq curated protein(タンパク質データ。TP53なら NP_000537 など)
    • NC_:Complete Genomic Molecule(ゲノム全体・染色体レベルのデータ)
    • XM_:Predicted mRNA(計算機で予測された未検証のmRNAデータ)
    つまり、NM_ とついているだけで「これは信頼性の高い、検証済みのmRNA(転写産物)データだな」と分かります。
  • 2. 連番 000546(識別ID):
    NM_ の後ろに続く6桁(または8桁)の数字は、特定の遺伝子・転写産物を一意に識別するための固有IDです。000546 は、ヒトの TP53遺伝子(主要な転写産物バリアント1) に割り当てられた専用の番号です。
  • 3. バージョン番号 .6(改訂履歴):
    ドット(.)の後ろの数字は、データの更新(リビジョン)回数を示しています。新しい実験データやゲノム解読技術の向上によって配列情報やアノテーション(注釈)が修正・更新されると、.1.2.3 … とカウントアップしていきます。.6 は、「TP53のこのmRNAデータが過去に5回修正され、現在第6版のデータである」ことを意味します。

② なぜTP53がサンプルに選ばれるのか?

この TP53 という遺伝子は、がん研究やバイオインフォマティクスにおいて「最も有名で、最も重要」と言っても過言ではない超重要遺伝子です。ITで言えば、「Hello World」 や 「サンプルコードで言うところの User クラス / main() 関数」 に相当する、バイオ界の超定番サンプルです。

  • 1. 「ゲノムの守護者(Guardian of the Genome)」と呼ばれる超有名遺伝子:
    TP53 は、細胞がガン化するのを防ぐがん抑制遺伝子です。DNAが傷ついた時に細胞周期を止めたり、修復不能な場合は細胞死(アポトーシス)を誘導したりする「エラー制御システム」の役割を果たします。
  • 2. ヒトのがんの半数以上で変異が見られる:
    あらゆる「がん研究」で真っ先に解析される遺伝子です。そのため、実験データ、論文、データベースの注釈(アノテーション)が世界一充実しており、チュートリアルや教科書で真っ先に取り上げられます。
  • 3. 配列長が手頃で学習に最適:
    転写産物(mRNA)の長さが 約 2,500 bp(塩基) と、短すぎず長すぎない非常に手頃なサイズです。プログラムのAPIで取得しても一瞬でダウンロードでき、メモリを圧迫せずにコードのテストができます。

2. NCBI APIからFASTAを取得するPythonコード

Biopython には、NCBIのWeb API(E-utilities)と通信するための Bio.Entrez モジュールが用意されています。これを使えば HTTP リクエストの構築を自作することなく、数行でデータを取り込めます。

from Bio import Entrez, SeqIO

# 1. NCBIマナーとしてメールアドレスを設定(必須)
Entrez.email = "your_email@example.com"

# 2. Entrez.efetch APIでNCBIからFASTAデータを取得
with Entrez.efetch(db="nuccore", id="NM_000546.6", rettype="fasta", retmode="text") as handle:
    # 3. レスポンスストリームを SeqIO.read で直にパース
    record = SeqIO.read(handle, "fasta")

# 4. 解析結果の確認
print(f"ID: {record.id}")
print(f"Description: {record.description}")
print(f"Sequence (先頭50文字): {record.seq[:50]}...")
print(f"Length: {len(record.seq)} bp")

3. 実行結果の確認

ID: NM_000546.6
Description: NM_000546.6 Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA
Sequence (先頭50文字): CTCAAAAGTCTAGAGCCACCGTCCAGGGAGCAGGTAGCTGCTGGGCTCCG...
Length: 2512 bp

4. コードとデータ処理のポイント

一連のコードの中で、Web APIとBiopythonがどのように連動しているかをエンジニア目線で解説します。

  • Entrez.email の設定規則:
    NCBIのAPIを利用する際は、利用者の識別用としてメールアドレスの設定が義務付けられています(過剰アクセス時の連絡用)。これを怠るとAPI制限や遮断の対象となるため、必ず設定しておきましょう。
  • Entrez.efetch()(HTTP GETリクエスト):
    NCBIのデータベース(db="nuccore")に対して、指定ID(NM_000546.6)のデータをFASTAテキスト形式(rettype="fasta", retmode="text")で問い合わせるAPI呼び出しです。
  • ストリームの直接パース:
    efetch が返すレスポンスハンドル(ストリーム)を、一旦ローカルファイルに保存することなくそのまま SeqIO.read(handle, "fasta") に流し込んでいます。これにより、メモリ上で完結するパイプラインが組めます。
  • 塩基長(2,512 bp):
    実際の実行結果から確認できる通り、正確に2,512 bpのmRNA配列が一瞬で取得され、SeqRecord オブジェクトとして扱えるようになります。スライス(record.seq[:50])を使って先頭配列(CTCAAAAG...)の確認も容易です。

5. ITエンジニア的まとめ

今回の実習のポイントは以下の通りです。

  • RefSeq IDの意味を理解する: NM_000546.6 のような文字列は、分子種(NM_)、固有ID(000546)、改訂版(.6)を示す意味のある識別子。
  • 公共DBからの自動取得: Bio.Entrez を使えば、手動ダウンロード不要でコード内から直接Web API経由でデータを取得できる。
  • ストリーム処理によるパイプライン化: APIレスポンスをそのまま SeqIO パーサーに渡して SeqRecord オブジェクト化可能。

これで「APIから実データを取得してプログラムに載せる」というデータインポートの手順が確立できました。

次回は【バイオインフォ実習】第5回:複数配列が含まれるマルチFASTAファイルのループ処理(SeqIO.parse)と統計解析に挑戦します!




【バイオインフォ実習】第3回:FASTA形式とBiopythonによる配列データ読み込み(SeqRecord)

前回はBiopythonを使ったセントラルドグマの基本操作を学びました。今回は、バイオインフォマティクスで最も頻繁に扱われるデータフォーマット「FASTA形式」と、それをBiopythonで読み込む基本コードを解説します。

配列データを扱うにあたり、まず「FASTA形式とは何か」をITエンジニア視点で整理しておきましょう。

1. FASTA形式とは?

FASTA(ファスタ)形式は、DNAやアミノ酸などの生物学的配列を扱うための標準的なテキストファイルフォーマットです。構造は至ってシンプルで、以下のように「ヘッダー」と「配列本体」で構成されています。

>seq1 Sample DNA sequence
ATGCGT

IT的に捉えると、以下のようなルールを持ったシンプルなデータ形式です。

  • ヘッダー行(メタデータ): 行頭が必ず > で始まります。> の直後の1単語が「配列ID」、それに続く文字列が「説明文(メタ情報)」となります。
  • 配列行(データ本体): 2行目以降にDNA塩基(A, T, G, C)が並びます。途中に改行が含まれていても、システム側では1つの連続したデータとして処理されます。
  • メタデータ付き構造化テキスト: JSONやXMLほど厳密ではありませんが、「ヘッダー+ペイロード」という構造を持った軽量なフォーマットです。

2. Biopythonでの実装コード

外部ファイルを用意する手間を省き、まずはコード内だけで自己完結する最小のFASTA読み込みプログラムを動かしてみましょう。io.StringIO を使うことで、文字列を疑似的なファイルとして読み込ませることができます。

from io import StringIO
from Bio import SeqIO

# 1. 最小構成のFASTA文字列を用意
fasta_data = """>seq1 Sample DNA sequence
ATGCGT"""


# 2. SeqIO.read でパース(単一配列の場合)
record = SeqIO.read(StringIO(fasta_data.strip()), "fasta")

# 3. 解析結果の出力
print(f"ID: {record.id}")
print(f"Description: {record.description}")
print(f"Sequence: {record.seq}")
print(f"Length: {len(record.seq)} bp")

3. 実行結果の確認

ID: seq1
Description: seq1 Sample DNA sequence
Sequence: ATGCGT
Length: 6 bp

4. record(SeqRecordクラス)とは何か?

コード内の record という変数には、一体どのような型が入っているのでしょうか?

このオブジェクトの正体は、Bio.SeqRecord.SeqRecord クラスのインスタンスです。Bio.SeqIO.read() 関数がパース処理を行った結果として返却されます。

  • 提示されている場所(定型パッケージ):
    Biopython 内の Bio.SeqRecord モジュールで定義されているクラスです。通常は直接インポートしなくても、Bio.SeqIO 経由でファイルを読み込んだ際に自動的に生成・返却されます。
  • 何のための型(オブジェクト)?:
    単なる塩基配列("ATGCGT" という文字列)だけでなく、「配列本体 +IDや説明文などのメタデータ」をセットで保持・管理するためのデータ構造(DTO/Entity)です。もし単なる文字列として配列を扱ってしまうと、「このDNAがどの遺伝子のものか」という重要情報が脱落してしまいます。それを防ぐため、配列(Seq オブジェクト)をコアに抱え込みつつ、属性情報(id, description など)を一つにまとめるカプセル化の役割を果たしています。

5. 出力結果の解説

プログラムを実行すると、FASTAテキストが自動的にパースされ、SeqRecord の各プロパティ(属性)として分解されたことがわかります。

  • ID: seq1record.id
    ヘッダー行(>seq1 Sample DNA sequence)から、> の直後にある最初の単語(スペース手前まで)が一意の識別子(ID)として自動抽出されます。ITで言えば主キー(Primary Key)のような扱いです。
  • Description: seq1 Sample DNA sequencerecord.description
    先頭の > を除いたヘッダー行全体のテキストが入ります。配列の名称や生物種情報などのメタデータがここに保持されます。
  • Sequence: ATGCGTrecord.seq
    改行などを除去した純粋な配列データ本体です。文字列のように見えますが、内部的には Biopython 独自の Seq オブジェクトになっており、相補鎖変換や転写・翻訳などのメソッドをそのまま呼び出せます。
  • Length: 6 bplen(record.seq)
    len() 関数を適用することで、配列の長さ(塩基数: base pairs)を簡単に取得できます。

6. ITエンジニア的まとめ

今回の実習のポイントは以下の通りです。

  • FASTAは標準のヘッダー付き配列フォーマット: > で始まるヘッダー+配列本体という、非常にシンプルなデータ構文。
  • SeqRecord はメタデータ保持コンテナ: 配列データだけでなく、IDや説明文を一つにまとめて管理する専用オブジェクト。
  • Bio.SeqIO は強力なパーサー: テキスト解析処理を自作することなく、オブジェクト(SeqRecord)として直感的にアクセスできる。

データフォーマットの構造さえ掴んでしまえば、バイオデータもファイルI/Oの1つに過ぎませんね。

次回は、複数の配列が含まれるFASTAファイルの読み込み(SeqIO.parse)や、NCBI等の公共データベースから実際のDNAデータを取得して解析する方法に挑戦します。


【生物学問題】老化と病気の原因を引き起こす!「活性酸素説」を攻略

好気性生物が代謝の過程で生成する反応性の高い酸素種が、生体高分子を酸化・損傷させ細胞老化を引き起こすとする説について整理しましょう。

1. 問題:活性酸素による細胞老化の仮説

【 問題 】 体内に存在する活性酸素(ROS)が、DNA・脂質・タンパク質などの生体高分子を酸化・損傷させることによって、細胞機能の低下や老化が進行するという考え方を何と呼ぶでしょうか?

① 活性酸素説   ② 限界寿命説   ③ 免疫低下説   ④ プログラム説

2. 正解:細胞老化のメカニズムに関する正解

正解: ① 活性酸素説(フリーラジカル説 / Free Radical Theory of Aging)

3. 解説:酸化ストレスとDNA損傷の蓄積

ミトコンドリアでの呼吸(ATP合成)の副産物として発生する活性酸素(ROS:Reactive Oxygen Species)は、強い酸化力を持ちます。

[ 活性酸素説のポイント ]
損傷対象:DNAの塩基改変(8-OHdGなどの酸化損傷)、タンパク質の変性、細胞膜脂質の過酸化脂質化。
修復限界と老化:生体にはSOD(スーパーオキシドジスムターゼ)やカタラーゼなどの抗酸化酵素系が存在しますが、修復しきれなかった損傷が長期間にわたって蓄積することで、細胞周期が恒久的に停止し(細胞老化)、機能不全に陥るとされます。
テロメア短縮との相互作用:ROSによるDNA二本鎖切断が、テロメア領域の短縮を加速させる要因にもなります。

1. 試験のポイント: 「体内に存在する活性酸素がDNAやたんぱく質を損傷させ、細胞を老化させる説=活性酸素説」です。ハーマン(Harmon)によって提唱された「フリーラジカル説」と同義である点も押さえておきましょう。
2. バイオインフォの視点: 酸化ストレス応答のトランスクリプトーム解析(RNA-Seq)では、ROS暴露下で発現が誘導されるストレス応答遺伝子群(NRF2シグナル経路など)の濃縮解析(GSEAやKEGG Pathway解析)が行われます。また、DNAシーケンスデータから酸化損傷由来の変異パターン(C>A/G>Tトランスバージョン等)のミューテーショナル・シグネチャーを定量解析する手法も用いられます。


4. まとめ

「活性酸素がDNAやたんぱく質を損傷させ、細胞を老化させる説=活性酸素説」です。システム開発で言えば、長期間の稼働(呼吸・代謝)に伴ってメモリリークやログファイルの肥大化(酸化損傷)が蓄積し、システム全体のパフォーマンスが徐々に低下して最終的にハングアップするような現象ですね!


【BI技術者認定試験対策】人間の染色体構成!「常染色体の数と構成」を攻略

ヒトの体細胞に含まれる染色体の種類とペア(対)の数について、常染色体と性染色体の組み合わせを含めて整理しましょう。

1. 問題:ヒトの染色体の種類と本数

【 問題 】 ヒト(体細胞)の常染色体は、[      ] 種類の染色体を [      ] 本ずつ持っている(計44本)。空欄に入る正しい数字の組み合わせはどれでしょうか?

① 22、2   ② 23、2   ③ 46、1   ④ 11、4

2. 正解:染色体数に関する正解

正解: ① 22、2

3. 解説:常染色体(22種×2本)と性染色体(1種×2本)

ヒトの体細胞には、合計で46本(23対)の染色体が存在します。これらは男女共通の「常染色体」と、性別を決定する「性染色体」に分かれます。

[ ヒトの染色体構成(合計46本 / 23対) ]
常染色体(Autosome):1番から22番までの22種類が存在し、父親由来と母親由来の相同染色体が2本ずつペアを形成(22種 × 2本 = 44本)。
性染色体(Sex chromosome):X染色体とY染色体の組合せ(女性:XX、男性:XY)で2本存在。
合計:常染色体44本 + 性染色体2本 = 46本(23対)

1. 試験のポイント: 問題文が「全体の染色体対(23種類を2本ずつ=46本)」を指しているのか、「常染色体のみ(22種類を2本ずつ=44本)」を指しているのかを冷静に読み取ることがポイントです。数字の22と23は引っかけとして頻出します。
2. バイオインフォの視点: ゲノムアライメント解析(BWAやSTARなど)では、リファレンスゲノムデータ(GRCh38など)として `chr1` ~ `chr22`(常染色体)、`chrX`、`chrY`(性染色体)、および `chrM`(ミトコンドリアDNA)のFASTAファイルを用いてマッピング処理が行われます。


4. まとめ

「ヒトは22種類の常染色体を2本ずつ持つ」です。システム開発で言えば、冗長化(ミラーリング)のために同じ規格のサーバー(染色体)を2台1組(22ペア)で運用し、さらにシステム種別判定用の特殊サーバー(性染色体)を2台追加してクラスタを組んでいるような構造ですね!