RNA-seq解析入門⑤：カウント編

この記事では、マッピング後のBAMファイルから遺伝子ごとのリード数を集計する「カウント」について解説します。BAMベースのfeatureCounts・HTSeq-countと、アラインメントフリーのSalmon・Kallistoの4ツールの特徴と使い分け、そしてHTSeq-countを使った実践コードの読み方まで説明します。

カウントとは？なぜ必要なのか

マッピングでリードをゲノム上に貼り合わせたら、次は「各遺伝子にいくつのリードが対応しているか」を数える作業が必要です。これをカウント（定量）と呼びます。

この数値が、後の発現変動解析（DESeq2・edgeRなど）で使う入力データになります。

💡 カウントとは？
ゲノム上の各遺伝子領域に重なっているリードの数を集計する作業です。この数値が大きいほど「その遺伝子がよく発現していた」ことを意味します。カウント結果は通常、遺伝子×サンプルの表形式（カウントマトリクス）として出力されます。

カウントの方法は大きく2つに分かれます。

📂 BAMベース：STARなどでマッピングしたBAMファイルを入力とする方法。featureCounts・HTSeq-countがこれにあたります。

⚡ アラインメントフリー：マッピングを行わず、FASTQファイルから直接カウントする方法。Salmon・Kallistoがこれにあたります。処理が非常に速いのが特徴です。

今回紹介するのパイプライン（STAR → HTSeq-count）はBAMベースの流れです。

FASTQファイル→トリミング→マッピング（STAR）→カウント（HTSeq-count）→発現変動解析

主なツールの紹介と比較

ツール名方式開発言語特徴
featureCountsBAMベースC高速・大規模サンプル向け
HTSeq-countBAMベースPythonシンプル・学習コストが低い
SalmonアラインメントフリーC++超高速・バイアス補正が充実
KallistoアラインメントフリーC++超高速・メモリ効率が良い

🗂 BAMベースのツール

featureCounts

SubreadパッケージのツールでC言語実装のため非常に高速です。数十サンプルを一括処理する際に特に力を発揮します。マルチスレッドにも対応しています。

メリット

  • 処理速度が非常に速い
  • 複数サンプルを一度にまとめて処理できる
  • マルチスレッド対応

デメリット

  • 細かいカウントモードの設定が必要な場面がある
  • HTSeqと比べてオプションが多く最初はやや難しい

このツールを選ぶ理由：サンプル数が多い解析や、パイプラインを自動化したい場合に最適です。

HTSeq-count

Python製のカウントツールで、コマンドがシンプルで動作が直感的に理解しやすいのが特徴です。RNA-seqの教育的なパイプラインでよく取り上げられています。

メリット

  • コマンドがシンプルで学習コストが低い
  • Python環境があればすぐ導入可能
  • 動作が透明で理解しやすい

デメリット

  • featureCountsより処理が遅い
  • 大規模サンプルには不向き

このツールを選ぶ理由：初学者やPythonに慣れている方に最適。解析の流れを丁寧に学びたい段階に向いています。

⚡ アラインメントフリーのツール

Salmon

FASTQファイルから直接トランスクリプトの発現量を推定するツールです。マッピングを省くため処理が超高速で、バイアス補正機能も充実しています。近年急速に普及しており、tximportと組み合わせてDESeq2に渡すワークフローが標準化されつつあります。

メリット

  • 処理速度が非常に速い（BAMベースの数十倍）
  • バイアス補正が充実している
  • トランスクリプトレベルの定量が可能

デメリット

  • BAMファイルが出力されないため中間ファイルの確認ができない
  • BAMベースと結果が異なる場合がある

このツールを選ぶ理由：処理速度と精度を両立させたい場合や、大規模サンプルをすばやく定量したい場面に向いています。

Kallisto

Salmonと同様にFASTQから直接定量するツールです。擬似アライメント（pseudo-alignment）という独自の手法を採用しており、非常に軽量・高速に動作します。

メリット

  • メモリ使用量が少ない
  • インストールが簡単
  • 処理速度が非常に速い

デメリット

  • Salmonと比べバイアス補正が限定的
  • BAMファイルが出力されない

このツールを選ぶ理由：メモリが限られた環境で高速に定量したい場合に向いています。

ツール選びのまとめ

状況おすすめの選択
はじめてカウントをするHTSeq-count
大量サンプルをBAMから速く処理したいfeatureCounts
マッピングを省いて高速に定量したいSalmon
メモリが少ない環境で高速定量したいKallisto

HTSeq-countの基本的な使い方

インストール

# pipでインストール pip install HTSeq # バージョン確認 htseq-count –version

bash 
# pipでインストール
pip install HTSeq

# バージョン確認
htseq-count --version

💡 Conda環境を使っている場合
conda install -c bioconda htseq でもインストールできます。

基本コマンド

STARで作成したソート済みBAMファイルとGTFアノテーションファイルを入力として使います。

bash 
htseq-count \
  -f bam \
  -r pos \
  -s no \
  -t exon \
  -i gene_id \
  output/sample_Aligned.sortedByCoord.out.bam \
  annotation.gtf \
  > counts.txt
-f bam入力ファイルの形式を指定します。BAMファイルを使う場合は bam を指定します。
-r posBAMファイルのソート方法を指定します。STARの出力はゲノム座標でソートされているため pos を指定します。
-s noストランド特異性の指定です。ストランド情報がない（非特異的）ライブラリの場合は no を指定します。ストランド特異的ライブラリでは yes または reverse に変更が必要です。
-t exonGTFファイルのどのfeature（特徴）をカウント対象にするかを指定します。通常はエクソン（exon）を指定します。
-i gene_idカウントをまとめる単位を指定します。gene_id を指定すると遺伝子単位でカウントが集計されます。
sample_Aligned…bamSTARが出力したBAMファイルを指定します。ファイル名はSTARの --outFileNamePrefix の設定によって変わります。
annotation.gtf遺伝子のアノテーション情報が入ったGTFファイルを指定します。マッピング時に使ったものと同じファイルを使うのが基本です。
> counts.txt結果を counts.txt というファイルに保存します。> はリダイレクト（出力先の指定）です。

💡 ストランド特異性（-s オプション）について
使用したライブラリ調製キットによって指定する値が変わります。不明な場合はRSeQCの infer_experiment.py で確認できます。一般的なTruSeq Stranded キットでは reverse を指定します。

出力結果のイメージ

実行後、counts.txt には以下のような形式で遺伝子ごとのカウント数が出力されます。

出力例（counts.txt） 
# 遺伝子ID              カウント数
ENSG00000000003    1523
ENSG00000000005       4
ENSG00000000419     892
ENSG00000000457     310
...
__no_feature       4821   # どの遺伝子にも対応しなかったリード
__ambiguous        1204   # 複数遺伝子に重なったリード
__too_low_aQual     302   # 品質スコアが低すぎたリード

💡 __no_feature などの特殊行について
カウント結果の末尾にはアンダースコアで始まる特殊な行が含まれます。これらはどの遺伝子にも割り当てられなかったリードの統計情報で、発現変動解析の前に除外する必要があります。

発展的なオプション（参考）

  • -m union複数の遺伝子に重なるリードの扱い方を指定（デフォルトは union）。intersection-strict や intersection-nonempty も選べる。
  • -a 10マッピング品質スコアの閾値を指定（デフォルト10）。低品質なマッピングを除外できる。
  • –additional-attr gene_name出力に遺伝子名（シンボル）列を追加する。結果が読みやすくなる。
  • -n 4使用するCPUコア数を指定。HTSeq v2.0以降で対応。
  • –counts-output counts.tsv> によるリダイレクトの代わりに出力ファイルを直接指定する方法。

まとめ

  • カウントは各遺伝子に対応するリード数を集計する作業で、発現変動解析の入力データになる
  • BAMベース（featureCounts・HTSeq-count）とアラインメントフリー（Salmon・Kallisto）の2方式がある
  • 初学者にはコマンドがシンプルなHTSeq-countがおすすめ
  • -s オプション（ストランド特異性）はライブラリキットに合わせて必ず確認する
  • 出力末尾の __no_feature などの特殊行は発現変動解析の前に除外する

次のステップ：カウントマトリクスが揃ったら、DESeq2またはedgeRを使って発現変動解析へと進みましょう。

コメント

タイトルとURLをコピーしました