Geminiのモデルの違い（RNA-seqの計算手法の質問）

結論

2.0 Flash = 2.0 Flash Thinking = Personalization  > Deep Research> 2.0 Pro 

2.0 Flash, Flash Thinking, Personalizationの違いがよくわからない。Deep Researchと2.0 Proはだめだった。


質問

FPKM, TPM, CPMをそれぞれ説明してください。さらにSTARでリファレンスゲノムにMappingした複数サンプルのBamファイルからそれぞれの値を計算し、サンプル間で比較するにはどのようなコマンドが必要か書いてください。

2.0 Flash

FPKM、TPM、CPMは、RNA-seqデータ解析において遺伝子発現量を定量化し、サンプル間で比較するために一般的に使用される正規化された測定値です。

各測定値の説明

  • FPKM (Fragments Per Kilobase of exon per Million mapped fragments):
    • マッピングされた断片数（ペアエンドシーケンスの場合はリードのペア数）を、遺伝子のエクソン長（キロベース）と全マッピング断片数（100万単位）で正規化した値です。
    • 遺伝子の長さとシーケンス深度を考慮しているため、異なる遺伝子間や異なるサンプル間で発現量を比較するのに適しています。
  • TPM (Transcripts Per Million):
    • まず遺伝子長でリードカウントを正規化し、次に全トランスクリプトの合計が100万になるように正規化した値です。
    • FPKMと同様に遺伝子の長さとシーケンス深度を考慮しますが、正規化の順序が異なります。
    • TPMは、サンプル間のトランスクリプト総量の違いによる影響を受けにくいため、サンプル間の遺伝子発現量比較に適しているとされています。
  • CPM (Counts Per Million mapped reads):
    • マッピングされたリード数を全マッピングリード数（100万単位）で正規化した値です。
    • シーケンス深度のみを考慮し、遺伝子の長さは考慮していません。
    • 主に、同じ遺伝子の発現量を異なるサンプル間で比較する場合に使用されます。

STARでマッピングしたBamファイルからの計算とサンプル間比較

STARでリファレンスゲノムにマッピングした複数サンプルのBamファイルからFPKM、TPM、CPMを計算し、サンプル間で比較するには、以下のツールとコマンドを使用できます。

  1. featureCounts (サブリード割り当て):
    • まず、Bamファイルから各遺伝子にマッピングされたリード数をカウントします。
Bash
featureCounts -p -t gene -g gene_id -a <GTFファイル> -o counts.txt <Bamファイル1> <Bamファイル2> ...
  • -p: ペアエンドリードを考慮します。
  • -t gene: 遺伝子レベルでカウントします。
  • -g gene_id: GTFファイル内の遺伝子ID属性を指定します。
  • -a <GTFファイル>: 遺伝子アノテーションのGTFファイルを指定します。
  • -o counts.txt: 出力ファイル名を指定します。
  • <Bamファイル1> <Bamファイル2> ...: 入力Bamファイルを指定します。
  1. RとDESeq2/tximport (FPKM/TPM/CPM計算とサンプル間比較):
    • featureCountsの出力ファイルをRに読み込み、DESeq2またはtximportパッケージを使用してFPKM、TPM、CPMを計算し、サンプル間比較を行います。

DESeq2を使用する場合のRコード例:

コード スニペット
# カウントデータを読み込む
counts <- read.table("counts.txt", header = TRUE, row.names = 1)

# サンプル情報を設定する
sampleTable <- data.frame(condition = factor(c("condition1", "condition2", ...)))

# DESeqDataSetを作成する
dds <- DESeqDataSetFromMatrix(countData = counts, colData = sampleTable, design = ~ condition)

# 正規化係数を計算する
dds <- estimateSizeFactors(dds)

# FPKM/TPMを計算する
fpkm <- fpkms(dds)
tpm <- tpms(dds)

# CPMを計算する
cpm <- counts(dds, normalized = TRUE)

# サンプル間比較 (DESeq2の機能を使用)
dds <- DESeq(dds)
results <- results(dds)

tximportを使用する場合のRコード例:

コード スニペット
# GTFファイルを読み込み、遺伝子長情報を取得する
library(rtracklayer)
gtf <- rtracklayer::import("<GTFファイル>")
txdb <- GenomicFeatures::makeTxDbFromGRanges(gtf)
tx2gene <- GenomicFeatures::select(txdb, keys = keys(transcripts(txdb)), columns = "GENEID", keytype = "TXID")

# featureCountsの出力ファイルを読み込む
files <- c("<counts.txt>")
names(files) <- c("sample1","sample2",...)
txi <- tximport::tximport(files, type = "gene", tx2gene = tx2gene)

# TPMを計算する
tpm <- txi$abundance

補足:

  • 上記のコマンドとRコードは一例であり、データや解析目的に応じて調整が必要です。
  • GTFファイルは、リファレンスゲノムに対応したアノテーションファイルを使用してください。
  • DESeq2は、統計的なサンプル間の比較に非常に有用なツールで、発現変動遺伝子の検出に利用可能です。
  • Rでの作業は、RとR studioをインストールすると容易に行えます。

これらの情報を参考に、RNA-seqデータの解析を進めてください。


2.0 Flash Thinking

はい、FPKM、TPM、CPMについてご説明し、STARでマッピングしたBAMファイルからこれらの値を計算してサンプル間で比較する方法について、必要なコマンド例と合わせてご説明します。

FPKM、TPM、CPMの説明

これらはRNA-Seqのリード数を正規化し、遺伝子発現量を定量化するための指標です。

  • FPKM (Fragments Per Kilobase of transcript per Million mapped reads): 主にシングルエンドRNA-Seqデータで使用される指標です。

    • 計算方法:
      1. 各遺伝子にマッピングされたリードの断片数（Fragments）をカウントします。
      2. そのカウント数を、遺伝子の長さ（キロベース）で割ります。これにより、遺伝子長による差が補正されます。
      3. さらに、全マッピングされたリードの総数（ミリオン）で割ります。これにより、ライブラリサイズによる差が補正されます。
    • 特徴: 転写産物の長さを考慮しているため、異なる長さの遺伝子間で発現量を比較するのに適しています。ただし、同じ遺伝子の異なるアイソフォーム間で発現量を比較するのには適していません。

  • TPM (Transcripts Per Million): シングルエンドおよびペアエンドRNA-Seqデータで使用できる指標です。近年ではFPKMよりも推奨されることが多いです。

    • 計算方法:
      1. 各遺伝子にマッピングされたリード（またはペアエンドの場合はフラグメント）数をカウントします。
      2. そのカウント数を、遺伝子の長さ（キロベース）で割ります。
      3. 上記で得られた値を全ての遺伝子について合計し、その合計値を100万で割ります。これにより、100万リードあたりの遺伝子長で正規化されたリード数が得られます。
      4. 最後に、各遺伝子の値を上記で得られた値で割ります。
    • 特徴: ライブラリサイズの影響を受けにくく、異なるサンプル間で遺伝子発現量を比較するのに適しています。また、遺伝子内の異なるアイソフォームの相対的な発現量をより正確に反映するとされています。

  • CPM (Counts Per Million mapped reads): シングルエンドおよびペアエンドRNA-Seqデータで使用できる指標です。遺伝子レベルでの発現量を比較するのに適しています。

    • 計算方法:
      1. 各遺伝子にマッピングされたリード（またはペアエンドの場合はフラグメント）数をカウントします。
      2. そのカウント数を、全マッピングされたリードの総数（ミリオン）で割ります。
    • 特徴: 遺伝子の長さを考慮していないため、異なる長さの遺伝子間で直接発現量を比較するのには適していません。主に、同じ遺伝子の異なるサンプル間での発現量の変化を見るのに適しています。また、非常に発現量の低い遺伝子に対して、FPKMやTPMよりも適切な場合があります。

BAMファイルから各値を計算し、サンプル間で比較する方法

STARでリファレンスゲノムにマッピングした後のBAMファイルからこれらの値を計算し、サンプル間で比較する一般的なワークフローは以下の通りです。

  1. リード数のカウント: 各遺伝子にマッピングされたリード数をカウントします。これには、featureCountsやHTSeq-countなどのツールがよく用いられます。
  2. 正規化と発現量の計算: カウントされたリード数をもとに、RのDESeq2やedgeR、tximportなどのパッケージを用いてFPKM、TPM、CPMを計算します。
  3. サンプル間の比較: 計算された発現量を用いて、サンプル間の遺伝子発現量の変化を解析します。DESeq2やedgeRなどのパッケージは、差分発現解析の機能も提供しています。

以下に、具体的なコマンド例を示します。ここでは、featureCountsでリード数をカウントし、RのDESeq2パッケージを用いて正規化と発現量の計算、そしてサンプル間の比較を行う例を示します。

前提条件:

  • STARによるマッピング後のBAMファイルが複数存在していること（例: sample1.bam, sample2.bam など）。
  • リファレンスゲノムのアノテーションファイル（GTFファイル）が存在すること。
  • Rがインストールされており、DESeq2パッケージがインストールされていること。

ステップ1: リード数のカウント

featureCountsを用いて各BAMファイルに対してリード数をカウントします。

Bash
# featureCountsのインストール (まだインストールしていない場合)
# 例: Linux環境の場合
sudo apt-get update
sudo apt-get install subread

# 各サンプルに対してリード数をカウント
featureCounts -a <path/to/your_annotation.gtf> -o sample1.counts.txt sample1.bam
featureCounts -a <path/to/your_annotation.gtf> -o sample2.counts.txt sample2.bam
# 他のサンプルも同様に実行
  • -a <path/to/your_annotation.gtf>: アノテーションファイルのパスを指定します。
  • -o <output_file>: 出力ファイル名を指定します。
  • <bam_file>: 入力となるBAMファイルを指定します。

複数のBAMファイルをまとめて処理することも可能です。

Bash
featureCounts -a <path/to/your_annotation.gtf> -o all_samples.counts.txt sample1.bam sample2.bam ...

この場合、出力ファイルは各サンプルごとのリード数を含むテーブル形式になります。

ステップ2: Rを用いた正規化と発現量の計算、サンプル間の比較

Rを起動し、DESeq2パッケージを用いて解析を行います。

コード スニペット
# DESeq2パッケージのインストール (まだインストールしていない場合)
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")

# ライブラリの読み込み
library(DESeq2)

# カウントデータの読み込み
# featureCountsで複数のBAMファイルをまとめて処理した場合
countData <- read.table("all_samples.counts.txt", header = TRUE, skip = 1, row.names = 1)
countData <- countData[, 6:ncol(countData)] # 必要な列（各サンプルのカウント数）を選択
colnames(countData) <- c("sample1", "sample2", ...) # 列名をサンプル名に変更

# サンプルのメタデータ（条件など）の作成
# 例: sample1とsample2がそれぞれcontrolとtreatmentの条件の場合
colData <- data.frame(condition = factor(c("control", "treatment", ...)))
rownames(colData) <- colnames(countData)

# DESeqDataSetの作成
dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = colData,
                              design = ~ condition)

# フィルタリング（低カウントの遺伝子を除去する場合）
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# 正規化と分散推定
dds <- DESeq(dds)

# FPKMとTPMの計算
fpkm_values <- fpkm(dds)
tpm_values <- tpm(dds)

# CPMの計算 (DESeq2には直接的な関数がないため、手動で計算)
cpm_values <- counts(dds, normalized = TRUE)

# 結果の確認 (例: sample1の最初の10遺伝子のFPKM, TPM, CPM)
head(fpkm_values[, "sample1"])
head(tpm_values[, "sample1"])
head(cpm_values[, "sample1"])

# サンプル間の比較 (例: treatment vs control)
results <- results(dds, contrast = c("condition", "treatment", "control"))
summary(results)

# 結果をテーブルとして保存
write.csv(as.data.frame(results), file="differential_expression_results.csv")
write.csv(fpkm_values, file="fpkm_values.csv")
write.csv(tpm_values, file="tpm_values.csv")
write.csv(cpm_values, file="cpm_values.csv")

ステップ3: 結果の比較

Rで計算されたFPKM、TPM、CPMの値が保存されたCSVファイルを読み込み、サンプル間で比較を行います。ヒートマップや散布図などを用いて可視化することも有効です。

コード スニペット
# 例: FPKM値を読み込んで最初の数行を表示
fpkm_data <- read.csv("fpkm_values.csv", row.names = 1)
head(fpkm_data)

# 例: TPM値をヒートマップで表示 (pheatmapパッケージが必要)
if (!requireNamespace("pheatmap", quietly = TRUE))
    install.packages("pheatmap")
library(pheatmap)
pheatmap(log2(tpm_data + 1), scale = "row") # log2変換して表示することが多い

# 例: 2つのサンプル間のTPM値を散布図で表示
plot(log2(tpm_data[, "sample1"] + 1), log2(tpm_data[, "sample2"] + 1),
     xlab = "log2(TPM + 1) of sample1", ylab = "log2(TPM + 1) of sample2")
abline(0, 1, col = "red") # y=xの線を追加

上記は基本的なワークフローの一例です。実際の解析では、実験デザインや目的に応じて適切なツールやパラメータを選択する必要があります。また、tximportパッケージを用いると、Salmonやkallistoなどのpseudo-alignmentツールで得られた定量結果をDESeq2に直接読み込むことも可能です。

ご自身の実験データや解析目的に合わせて、上記コマンドやRスクリプトを適宜修正してご活用ください。


Personalization

FPKM、TPM、CPMについてそれぞれ説明します。また、STARでリファレンスゲノムにマッピングした複数サンプルのBamファイルからそれぞれの値を計算し、サンプル間で比較するための一般的なワークフローと関連ツールについて説明します。

FPKM (Fragments Per Kilobase of transcript per Million mapped reads)

  • 遺伝子発現量を定量化するための指標の一つで、主にペアエンドRNA-seqデータで使用されます。
  • 100万個のマッピングされたリードあたりの、遺伝子の長さ1kbあたりのリードペア数として定義されます。
  • 遺伝子の長さとライブラリサイズ（総リード数）を正規化することで、異なる長さの遺伝子や異なるシーケンス深度のサンプル間での比較をある程度可能にします。

TPM (Transcripts Per Million)

  • 近年、FPKMよりも推奨されることの多い正規化手法で、シングルエンドとペアエンドの両方のRNA-seqデータで使用できます。
  • 以下の2段階で計算されます。
    1. まず、1000個のヌクレオチドあたりのリード数を各遺伝子について計算します。これは、各遺伝子のリード数をその遺伝子の長さに1000を掛けたもので割ることで得られます。
    2. 次に、ステップ1で得られた値を全ての遺伝子で合計し、その合計値で各遺伝子の値を割って100万を掛けます。
  • TPMの利点は、全ての遺伝子のTPM値の合計が100万になるため、サンプル間の比較がより直感的になることです。

CPM (Counts Per Million mapped reads)

  • シングルエンドRNA-seqデータで一般的に使用される正規化手法ですが、ペアエンドデータにも適用できます。
  • 100万個のマッピングされたリードあたりのリード数として定義されます。
  • 遺伝子の長さによる正規化は行われないため、異なる長さの遺伝子間で直接比較することは適切ではありませんが、同じ遺伝子の異なるサンプル間での発現量の変化を見るのに適しています。

STARでマッピングした複数サンプルのBamファイルからそれぞれの値を計算し、サンプル間で比較する方法

STARはリードをゲノムやトランスクリプトームにアラインするツールであり、FPKM、TPM、CPMの値を直接出力するわけではありません。これらの値を計算するには、STARのアウトプット（通常はBAMファイル）を別のツールに入力する必要があります。

一般的なワークフローは以下の通りです。

  1. リードカウントの取得:

    • STARのアラインメント後、各遺伝子にマッピングされたリード数をカウントする必要があります。これには、featureCounts（Rsubreadパッケージ）やhtseq-count（HTSeqパッケージ）などのツールがよく用いられます。これらのツールは、BAMファイルを遺伝子のアノテーション（通常はGTFファイル）と照らし合わせ、各遺伝子にマッピングされたリード数をカウントします。
    Bash
    # featureCountsの例
featureCounts -a annotation.gtf -t exon -g gene_id -o counts.txt sample1.bam sample2.bam ...

# htseq-countの例
htseq-count -f bam sample1.bam annotation.gtf
htseq-count -f bam sample2.bam annotation.gtf

  2. FPKM、TPM、CPMの計算とサンプル間の比較:

    • リードカウントが得られたら、Rなどの統計解析環境でこれらの値を計算し、サンプル間で比較します。よく使われるRのパッケージとしては、DESeq2やedgeRなどがあります。これらのパッケージは、リードカウントデータを受け取り、正規化（FPKM、TPM、CPM相当の値を計算可能）、差分発現解析などを実行できます。

    以下は、RとDESeq2を用いた一般的な手順の例です。

    コード スニペット
    # DESeq2パッケージのインストールとロード
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")
library(DESeq2)

# カウントデータの読み込み (featureCountsなどの出力ファイルを読み込む)
countData <- read.table("counts.txt", header=TRUE, row.names=1)

# サンプルの情報を作成
sampleInfo <- data.frame(condition=factor(c("control", "control", "treatment", "treatment"))) # 例

# DESeqDataSetオブジェクトの作成
dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = sampleInfo,
                              design = ~ condition)

# 正規化と分散の推定
dds <- DESeq(dds)

# 正規化されたカウントの取得 (TPMに類似したスケールされた値)
normalized_counts <- counts(dds, normalized=TRUE)
head(normalized_counts)

# FPKMまたはTPMの計算 (DESeq2のデフォルト出力には直接含まれていませんが、以下の方法で計算できます)
# 注: 遺伝子の長さ情報が必要です (GTFファイルから取得するなど)。
gene_lengths <- ... # 各遺伝子の長さを格納したベクトル

# FPKMの計算 (DESeq2のnormalized_countsを使用)
fpkm <- apply(normalized_counts, 2, function(x) (x * 1e9) / (gene_lengths / 1000) / sum(x))
head(fpkm)

# TPMの計算 (DESeq2のnormalized_countsを使用)
tpm <- apply(countData, 2, function(x) {
  rate <- x / (gene_lengths / 1000)
  rate / sum(rate) * 1e6
})
head(tpm)

# CPMの計算 (DESeq2のnormalized_countsを使用)
cpm <- apply(countData, 2, function(x) (x / sum(x)) * 1e6)
head(cpm)

# サンプル間の比較や差分発現解析
results <- results(dds)
summary(results)

補足事項:

  • 遺伝子の長さの情報は、使用しているゲノムアノテーションファイル（GTFファイルなど）から取得する必要があります。
  • サンプル間の比較を行う際には、適切な正規化手法を選択することが重要です。DESeq2やedgeRは、より高度な正規化手法を提供しており、単純なFPKM、TPM、CPMだけでなく、統計的なモデルに基づいた差分発現解析を行うことができます。
  • 大規模なデータセットの場合、コマンドラインツールとRなどの統計環境を組み合わせることで、効率的な解析が可能です。
  • 最近では、Salmonやkallistoといったpseudo-alignment/quantificationツールも広く使われています。これらのツールはアラインメントと定量化をより高速に行うことができ、TPMなどの値を直接出力します。

ユーザーの検索履歴から、screpertoireや免疫関連のキーワードが見られるため、シングルセルRNA-seqデータに関心がある可能性もあります。シングルセルRNA-seqの場合、正規化や比較の手法はバルクRNA-seqとは異なる場合があります。しかし、基本的なリードカウントからFPKM、TPM、CPM相当の値を計算する概念は共通しています。

具体的なコマンドは、使用するツールやデータの形式、解析の目的に応じて変わってきますので、上記は一般的なガイドラインとして参考にしてください。必要であれば、具体的なツールや解析についてより詳細な情報を提供できます。

2.0 Pro

はい、承知いたしました。FPKM, TPM, CPM の説明と、STARでマッピングしたBAMファイルからの計算方法、サンプル間比較について解説します。

1. FPKM, TPM, CPM の定義

これらはすべて、RNA-Seq データにおいて、遺伝子発現量を定量化するための正規化手法です。それぞれの指標は、異なる方法でリードカウントを補正し、遺伝子長やサンプル間のシーケンス深度の違いによる影響を軽減します。

  • FPKM (Fragments Per Kilobase of exon per Million mapped fragments):

    • 定義: 100万個のマッピングされたフラグメントあたり、かつ、対象遺伝子のエクソン領域1キロベースあたりのフラグメント数。
    • 特徴:
      • ペアエンドRNA-Seqデータで用いられることが多い（フラグメント単位でカウントするため）。
      • 遺伝子長とライブラリサイズ（総リード数）の両方で補正される。
      • 同一サンプル内での遺伝子間比較に適している。
      • 異なるサンプル間での比較には、ライブラリサイズの違いが完全に補正されないため、適さない場合がある。

  • TPM (Transcripts Per Million):

    • 定義: 100万個の転写産物あたり、対象遺伝子の転写産物数。
    • 特徴:
      • 遺伝子長で補正した後、サンプル内の全転写産物の総数で割る。
      • TPMの合計はサンプル間で一定（100万）になるため、サンプル間比較に適している。
      • 遺伝子長の長い遺伝子ほど、同じ発現量でもリードカウントが多くなるバイアスを考慮できる。
      • 比較的発現量の多い遺伝子に有利になる傾向がある。

  • CPM (Counts Per Million):

    • 定義: 100万個のマッピングされたリードあたり、対象遺伝子にマッピングされたリード数。
    • 特徴:
      • 最も単純な正規化手法。ライブラリサイズ（総リード数）のみで補正。
      • 遺伝子長は考慮されない。
      • サンプル間での比較は、遺伝子長が同じ遺伝子に限定する場合や、遺伝子長の影響を無視できる場合に適している。
      • 発現変動の大きな遺伝子を検出するのに適している場合がある。

2. STARでマッピングしたBAMファイルからの各値の計算

STAR は、RNA-Seq リードをリファレンスゲノムにマッピングするためのアライナーです。STAR の出力である BAM ファイルには、各リードがゲノムのどの位置にマッピングされたかの情報が含まれています。

各値 (FPKM, TPM, CPM) を計算するには、一般的に以下の手順を踏みます。

  1. リードカウント: 各遺伝子（または転写産物）にマッピングされたリード数をカウントします。
  2. 正規化: カウントしたリード数を、遺伝子長やライブラリサイズで補正し、FPKM, TPM, または CPM を算出します。

これらの計算を自動化するツールがいくつかあります。

  • RSEM: 転写産物レベルでの発現定量に特化したツール。TPM と FPKM を計算できます。

    • STARと連携して使用することが一般的です。
    Bash
    # STAR でマッピング (例)
STAR --runThreadN 12 --genomeDir /path/to/genome/index --readFilesIn read1.fastq.gz read2.fastq.gz --outSAMtype BAM SortedByCoordinate --outFileNamePrefix sample1_

# RSEM で定量 (例)
rsem-calculate-expression --bam --paired-end sample1_Aligned.sortedByCoord.out.bam /path/to/rsem/reference sample1

  • featureCounts (Subread package): 遺伝子レベルでのリードカウントを行うツール。

    • 高速で、大量のBAMファイルを処理できます。
    • カウント結果を元に、Rなどの統計解析ツールでFPKM, TPM, CPMを計算できます。
    Bash
    featureCounts -a genes.gtf -o counts.txt -p -t exon -g gene_id sample1.bam sample2.bam sample3.bam

  • StringTie: 転写産物のde novoアセンブリと発現定量を行うツール

    • アセンブリと同時にTPM、FPKMの計算が可能。
    Bash
    stringtie sample1.bam -p 12 -G genes.gtf -o sample1.gtf -l sample1 -A gene_abundances_sample1.tsv

  • Salmon: 準マッピング(alignment-free)に基づく高速な発現定量ツール。TPM を計算できます。

    Bash
    salmon quant -i /path/to/salmon/index -l A -1 read1.fastq.gz -2 read2.fastq.gz -o sample1_quant -p 12

3. サンプル間比較

計算された FPKM, TPM, または CPM 値を用いて、サンプル間で遺伝子発現量を比較します。一般的には、以下の方法が用いられます。

  • R などの統計解析ソフトウェア:

    • DESeq2, edgeR などのパッケージを利用して、正規化と統計的検定 (例: t検定, 分散分析) を行い、有意な発現変動を示す遺伝子を検出します。
    • これらのパッケージは、CPMや生のカウントデータを入力として、独自の正規化（TMM正規化など）を行い、より精密な比較を可能にします。
    • ヒートマップや火山プロットなどの可視化も容易です。
    コード スニペット
    # DESeq2 の例 (カウントデータから)
library(DESeq2)
# ... (カウントデータの読み込み、実験デザインの設定など) ...
dds <- DESeqDataSetFromMatrix(countData = counts, colData = sample_info, design = ~ condition)
dds <- DESeq(dds)
results <- results(dds)

  • 専用のツール:

    • sleuth は、kallisto や Salmon の出力 (TPM) を直接入力として、サンプル間比較を行うことができます。

重要な注意点

  • 解析の目的に応じて適切な指標を選ぶ:
    • サンプル内での遺伝子間比較: FPKM または TPM
    • サンプル間での比較: TPM または CPM (遺伝子長が同じ場合)
    • 発現変動の大きい遺伝子を検出: CPM
    • より精密な比較: DESeq2 や edgeR などで、生のカウントデータから正規化と検定を行う
  • 遺伝子長の情報: FPKM と TPM の計算には、遺伝子長 (またはエクソン領域の長さ) の情報が必要です。GTF/GFF ファイルから取得できます。
  • マルチマッピングリードの扱い: 複数の場所にマッピングされるリード (マルチマッピングリード) の扱いは、ツールによって異なります。
  • 実験デザイン: サンプル間比較を行う際は、実験デザイン（例：対照群と処理群）を考慮し、適切な統計手法を選択することが重要です。

上記を参考に、目的に合ったツールと指標を選択し、RNA-Seq データの解析を進めてください。

Deep Research

RNA-SeqにおけるFPKM、TPM、CPMの解説とサンプル間比較のためのコマンド

1. RNA-Seqと遺伝子発現定量化の概要

RNA-Seq（RNAシーケンス）は、次世代シーケンサー（NGS）技術を用いて遺伝子発現を検出する手法です。この技術は、RNA抽出から始まり、mRNAの断片化、cDNAへの逆転写、そしてシーケンスライブラリの作成を経て、数百万ものリード（短いDNA配列）を生成します。これらのリードは、参照ゲノムまたはトランスクリプトームにアラインメントされ、cDNAライブラリの質的および量的な組成を明らかにします 1。RNA-Seqの主な目的の一つは、RNA転写産物の存在量を測定し、それによって遺伝子活動を評価することです。初期段階で得られるローリードカウントは、遺伝子の発現レベルを示す指標となりますが、サンプルの調製過程で生じる技術的なばらつき（RNA量、品質、増幅効率の違いなど）の影響を受けやすく、そのままでは生物学的な差異を正確に反映しているとは言えません 2.

遺伝子発現の正規化は、これらの技術的なばらつきを調整し、真の生物学的な違いをより正確に捉えるために不可欠なプロセスです 2。正規化が必要となる主な要因としては、シーケンス深度（サンプルあたりの総リード数の違い）4、遺伝子長のばらつき（長い遺伝子ほど多くのリードがマッピングされやすい）4、そしてRNA組成の偏り（一部の高度に発現する遺伝子が他の遺伝子の相対的な発現量に影響を与える）5 などが挙げられます。正規化の目的は、カウントマトリックスを相対的な値に変換し、技術的なバイアスを取り除き、その後の統計的な差異解析をより有意なものにすることにあります 2。

本レポートでは、RNA-Seqデータ解析において一般的に用いられる正規化手法であるFPKM（Fragments Per Kilobase Million）、TPM（Transcripts Per Million）、およびCPM（Counts Per Million）について詳細に解説します。それぞれの定義、計算方法、生物学的な解釈、そして利用場面について説明するとともに、STARアライメントによって得られたBAMファイルからこれらの値を算出し、複数のサンプル間で比較するための具体的なコマンドライン操作についても解説します。RNA-Seqデータの正規化手法の進化は、データ解析の精度と解釈性を向上させるための継続的な取り組みを示しており、初期の方法が基本的なバイアスに対処したのに対し、より新しい方法は組成バイアスのようなより複雑な問題にも対応しています。

2. FPKM（Fragments Per Kilobase Million）の詳細解説

2.1. 定義と生物学的解釈

FPKMは、Fragments Per Kilobase of exon model per Million mapped fragmentsの略であり、主にペアエンドRNA-Seqデータのために開発された正規化手法です 2。ペアエンドシーケンスでは、一つのDNA断片の両端からリードが生成されるため、「フラグメント」という単位が用いられます。シングルエンドRNA-Seqデータに対する同様の指標としては、RPKM（Reads Per Kilobase Million）があります 1。FPKMは、遺伝子長とシーケンス深度の両方を考慮して正規化された遺伝子の相対的な発現レベルを示します。高いFPKM値は、その遺伝子の発現量が多いことを意味します 2。

2.2. 段階的な計算と数式

FPKMの計算は、以下の数式に基づいて行われます 9。

FPKM = (特定の遺伝子にマッピングされたフラグメント数 × 10<sup>9</sup>) / (全マッピングされたリード数 × 遺伝子長（塩基対）)

この計算は、以下の2つの主要な正規化ステップを含んでいます。

  1. 遺伝子長の正規化: まず、特定の遺伝子にマッピングされたフラグメント数をその遺伝子の長さ（キロベース）で割ります。これは、長い遺伝子ほど多くのフラグメントが生成される可能性が高いという事実を考慮するためです 2。
  2. シーケンス深度の正規化: 次に、長さで正規化されたカウントを、全マッピングされたフラグメント数（ミリオン単位）で割ります。これにより、ライブラリサイズの差異が補正されます 2。

2.3. 利用場面と限界

FPKMは、同一サンプル内の異なる遺伝子間の発現量を比較したり、異なるサンプル間で特定の遺伝子の発現量を比較したりする際に用いられます 2。技術的なレプリケート間の比較にも有用です 8。

しかしながら、FPKMにはいくつかの限界も存在します。主なものとして、RNA組成バイアスが挙げられます 3。FPKMは、全RNA出力がすべてのサンプルで同じであると仮定していますが、これは必ずしも真実ではありません。一部の高度に発現する遺伝子が全リード数に大きな影響を与え、他の遺伝子の発現量の解釈を誤らせる可能性があります 3。また、多くの差分発現解析ツールは、正規化されたカウントではなくローカウントを必要とするため、FPKMは差分発現解析には理想的ではありません 8。さらに、サンプル内の全遺伝子のFPKM値の合計は一定ではないため、異なるサンプル間で発現量の割合を直接比較することが難しいという点も課題です 10。FPKMにおける正規化の順序（深度を最初に、次に長さを考慮）が、サンプル間の比較における限界の一因となっています。全リード数で最初に正規化することで、その後の長さの正規化が各サンプルのトランスクリプトーム全体の組成に依存するようになります。もしあるサンプルに非常に高度に発現する遺伝子がいくつか存在する場合、シーケンス深度の正規化係数が大きくなり、そのような高発現遺伝子を持たないサンプルと比較して、他の遺伝子のFPKM値が相対的に低くなる可能性があります。

3. TPM（Transcripts Per Million）の詳細解説

3.1. 定義と生物学的解釈

TPMは、Transcripts Per Kilobase of exon model per Million mapped readsの略であり 2、転写産物の相対的な存在量を示す指標です。サンプル内のすべての転写産物のTPM値の合計は常に100万になるため、各転写産物がトランスクリプトーム全体に占める割合として解釈することができます 2。

3.2. 段階的な計算と数式

TPMは、以下のいずれかの方法で計算できます。

方法1:

  1. RPK（Reads Per Kilobase）の計算: 各遺伝子について、マッピングされたリード数をその遺伝子の長さ（キロベース）で割ります 6。 RPK = 特定の遺伝子にマッピングされたリード数 / 遺伝子長（キロベース）
  2. 全RPK値の合計: サンプル内のすべての遺伝子のRPK値を合計します。
  3. TPMの計算: 各遺伝子のRPK値を全RPK値の合計で割り、100万を掛けます 6。 TPM = (特定の遺伝子のRPK / 全遺伝子のRPKの合計) × 1,000,000

方法2（FPKMから変換）:

  1. サンプル内の全遺伝子のFPKM値を合計します。
  2. 各遺伝子のFPKM値を全FPKM値の合計で割り、100万を掛けます 14。 TPM<sub>i</sub> = (FPKM<sub>i</sub> / Σ<sub>j</sub> FPKM<sub>j</sub>) × 1,000,000

TPMの計算では、最初に遺伝子長で正規化し 2、次いでサンプル内の長さ正規化されたカウントの合計で正規化します 2。

3.3. 利用場面とFPKMに対する利点

TPMは、サンプル内およびサンプル間の遺伝子発現レベルの比較に適しています 2。

FPKMに対する主な利点は、異なるサンプル間でより正確な比較ができることです 2。TPMでは、サンプル内の全TPM値の合計が常に100万であるため、各転写産物の割合を直接比較できます。また、TPMは、最初に遺伝子長で正規化することで、高度に発現する遺伝子が他の遺伝子の正規化された値に与える影響を軽減します 3。FPKMの開発者自身も、現在ではTPMを推奨しています 12。TPMにおける正規化の順序（長さを最初に、次に深度を考慮）が、FPKMと比較してサンプル間の比較性能が向上する主要な理由です。遺伝子長で最初に正規化することで、各サンプルのトランスクリプトーム内における各転写産物の相対的な存在量が、全体のシーケンス深度や他の遺伝子の発現レベルに依存することなく、より正確に表現されるようになります。これにより、異なるサンプル間で正規化された値をより適切に比較することが可能になります。

4. CPM（Counts Per Million）の詳細解説

4.1. 定義と生物学的解釈

CPMは、Counts Per Million mapped readsの略であり 1、サンプル内の全リード数（ミリオン単位）で正規化された、特定の遺伝子にマッピングされたリード数を示します。主にシーケンス深度の違いを補正するために用いられます 1。

4.2. 段階的な計算と数式

CPMは、以下の数式で計算されます 1。

CPM = (特定の遺伝子にマッピングされたリード数 / サンプル内の全マッピングされたリード数) × 1,000,000

この計算では、各遺伝子のローカウントをサンプル内の全マッピングされたリード数で割り、その結果に100万を掛けることで、シーケンス深度のみを考慮した正規化を行います 1。

4.3. 利用場面と限界

CPMは、シーケンス深度の単純な正規化として利用されます 1。また、低発現遺伝子をフィルタリングする際（CPMの閾値を使用）8 や、特定の種類の可視化にも用いられます 7。

しかしながら、CPMは遺伝子長を考慮しないため 3、同一サンプル内の異なる遺伝子間の発現量を比較するのには適していません 3。また、FPKMと同様に、RNA組成バイアスにも影響を受けやすいです 3。CPM単独では、差分発現解析には推奨されません。より高度な手法が、データの分布やばらつきを考慮するために用いられます 8。CPMは、議論されている中で最も基本的な正規化手法であり、主にライブラリサイズの差異に対処します。遺伝子長の正規化を行わないため、多くの下流解析での有用性は限定的です。長い遺伝子は本質的により多くのリードを引き付けるため、CPM値は異なる長さの遺伝子を比較する際に誤解を招く可能性があります。例えば、短く高度に発現する遺伝子が、長く中程度に発現する遺伝子と類似したCPMを持つことがあります。

5. STARアラインメントされたBAMファイルからのFPKM、TPM、CPMの計算

5.1. 前提条件

  • 各サンプルのSTARアラインメントされたBAMファイル。
  • 遺伝子座標とエキソン情報を含むGTFまたはGFF形式の遺伝子アノテーションファイル。

5.2. 段階的なコマンドラインアプローチ

以下の手順では、一般的なコマンドラインツールを用いてFPKM、TPM、およびCPMを計算する方法を示します。

ステップ1: featureCountsを用いたロー遺伝子カウントの取得

featureCountsは、Subreadパッケージに含まれるツールで、ゲノム特徴（ここでは遺伝子）にマッピングされたリード数をカウントするために広く使用されています。

Bash
featureCounts -a annotation.gtf -t exon -g gene_id -o gene_counts.txt sample1.bam sample2.bam ...
  • -a annotation.gtf: アノテーションファイルを指定します。
  • -t exon: カウントするフィーチャータイプとしてエキソンを指定します。
  • -g gene_id: GTFファイル内の遺伝子を定義する属性を指定します。
  • -o gene_counts.txt: 遺伝子カウントの出力ファイル名を指定します。
  • sample1.bam sample2.bam ...: 入力となるBAMファイルを指定します。

出力ファイル gene_counts.txt には、遺伝子IDと各サンプルにおけるローリードカウントを含むテーブルが生成されます。

ステップ2: アノテーションファイルからの遺伝子長の抽出

FPKMおよびTPMの計算には遺伝子長が必要です。遺伝子長は通常、その遺伝子に属するすべてのエクソンの長さの合計として定義されます。以下の awk コマンドを使用して、GTFファイルからエクソン座標を抽出し、遺伝子長を計算できます。

Bash
awk 'BEGIN{FS="\t"} $3=="exon" {lengths[$9] += $5 - $4 + 1} END{for (gene in lengths) {gsub(/gene_id \"|\";/, "", gene); print gene "\t" lengths[gene]/1000}}' annotation.gtf > gene_lengths.txt

このスクリプトは、GTFファイルから各遺伝子のエクソン座標を抽出し、それらの長さを合計し、キロベース単位に変換します。出力は、遺伝子IDとその長さを含むタブ区切りファイル gene_lengths.txt です。

ステップ3: awkを用いたCPMの計算

CPMの計算に必要なのは、ローカウントと総マッピングリード数のみです。まず、各BAMファイルの総マッピングリード数を計算します（この情報が featureCounts の出力に直接含まれていない場合や、クロスチェックとして）。

Bash
samtools view -c sample1.bam
samtools view -c sample2.bam
...

次に、以下の awk コマンドを使用してCPMを計算します。

Bash
awk 'BEGIN{FS="\t"; OFS="\t"} NR==1 {print $0, "CPM_"FILENAME} NR>1 {sum=0; for(i=2; i<=NF; i++) {sum+=$i}; for(i=2; i<=NF; i++) {cpm = ($i / sum) * 1000000; printf "%s%s%.2f%s", $i, OFS, cpm, (i==NF ? "\n" : OFS)}}' gene_counts.txt

このスクリプトは、遺伝子カウントファイルを読み込み、各サンプルの総リード数を計算し、各遺伝子のCPM値を算出します。

ステップ4: awkを用いたFPKMの計算

まず、遺伝子カウントと遺伝子長を結合します。

Bash
join -t $'\t' gene_counts.txt gene_lengths.txt > merged_counts_lengths.txt

次に、以下の awk コマンドを使用してFPKMを計算します。

Bash
awk 'BEGIN{FS="\t"; OFS="\t"} NR==1 {print $0, "FPKM_"FILENAME} NR>1 {length=$NF; total_reads[$1]=0; for(i=2; i<NF; i++) {total_reads[$1]+=$i}; for(i=2; i<NF; i++) {fpkm = ($i / (length * (total_reads[$1]/1000000))); printf "%s%s%.2f%s", $i, OFS, fpkm, (i==NF-1 ? "\n" : OFS)}}' merged_counts_lengths.txt

このスクリプトは、結合されたファイルから遺伝子長とカウントを読み取り、各遺伝子のFPKM値を計算します。

ステップ5: awkを用いたTPMの計算

merged_counts_lengths.txt ファイルを使用して、まずRPKを計算します。

Bash
awk 'BEGIN{FS="\t"; OFS="\t"} NR==1 {print $0, "RPK_"FILENAME} NR>1 {length=$NF; for(i=2; i<NF; i++) {rpk[i] = $i / length; printf "%s%s%.2f%s", $i, OFS, rpk[i], (i==NF-1 ? "\n" : OFS)}}' merged_counts_lengths.txt > gene_rpk.txt

次に、以下のコマンドでTPMを計算します。

Bash
awk 'BEGIN{FS="\t"; OFS="\t"} NR==1 {print $0, "TPM_"FILENAME} NR>1 {total_rpk=0; for(i=2; i<=NF; i++) {total_rpk+=$i}; for(i=2; i<=NF; i++) {tpm = ($i / total_rpk) * 1000000; printf "%s%s%.2f%s", $i, OFS, tpm, (i==NF ? "\n" : OFS)}}' gene_rpk.txt

この2段階のプロセスにより、各遺伝子のTPM値が得られます。

表1: ローカウント、遺伝子長、および計算された正規化値の例

遺伝子ID遺伝子長 (kb)サンプル1 ローカウントサンプル2 ローカウントサンプル1 CPMサンプル2 CPMサンプル1 FPKMサンプル2 FPKMサンプル1 TPMサンプル2 TPM
geneA1.515003000150.00300.00100.00200.00130.00260.00
geneB3.030006000300.00600.00100.00200.00130.00260.00
geneC0.75750150075.00150.00100.00200.00130.00260.00
geneD2.020004000200.00400.00100.00200.00130.00260.00

上記の表は、架空のデータを用いた正規化値の例を示しています。実際のデータでは、より多くの遺伝子とサンプルが含まれることになります。コマンドラインによるアプローチは、正規化の基礎となる計算を理解する上で重要ですが、大規模なデータセットに対しては、専用のRパッケージを使用する方が効率的で堅牢な解析につながります。

6. 複数のサンプル間でのFPKM、TPM、CPM値の比較

6.1. 計算された値の整理

前のステップで得られたタブ区切りファイルは、Rのような統計ソフトウェアに簡単にインポートしたり、コマンドラインツールでさらに処理したりできます。

6.2. コマンドライン戦略による値の比較

  • awkとgrepを用いた基本的な比較とフィルタリング:

    • 特定の遺伝子のTPM値をすべてのサンプルで検索する場合：
      Bash
      grep "gene_name" gene_tpm.txt
    • サンプル内の全遺伝子の平均CPMを計算する場合：
      Bash
      awk 'BEGIN{FS="\t"; sum=0; count=0} NR>1 {sum += $2; count++} END{if(count>0) print "Average CPM:", sum/count}' gene_cpm.txt

  • Rを用いた統計解析の導入:

    • Rは、RNA-Seqデータの統計解析と可視化のための強力なツールを提供します。差分発現解析には通常、edgeR 1 や DESeq2 8 のようなパッケージを用いてローカウントが用いられますが、正規化された値（特にTPMとCPM）は、全体的な傾向を可視化し探索するのに役立ちます。

  • 基本的なRコマンドによるデータの読み込みと比較:

    • CPMデータをRに読み込む：
      コード スニペット
      cpm_data <- read.table("gene_cpm.txt", header=TRUE, sep="\t")
    • 特定の遺伝子の平均CPMをサンプル間で計算する：
      コード スニペット
      gene_name <- "your_gene_id"
mean_cpm <- mean(as.numeric(cpm_data[cpm_data$gene_id == gene_name, 2:ncol(cpm_data)]))
print(mean_cpm)
    • TPM値のサンプル間の分布を箱ひげ図で表示する：
      コード スニペット
      tpm_data <- read.table("gene_tpm.txt", header=TRUE, sep="\t")
boxplot(tpm_data[, 2:ncol(tpm_data)], main="TPM Distribution Across Samples", ylab="TPM")
    • TPM値に基づいてサンプル間の相関行列を計算する：
      コード スニペット
      cor_matrix <- cor(t(tpm_data[, 2:ncol(tpm_data)]))
print(cor_matrix)

コマンドラインツールは基本的なデータ操作には役立ちますが、複数のサンプルにわたるRNA-Seqデータを解釈するためには、統計解析と可視化のためのより包括的な環境を提供するRが不可欠です。サンプル間の遺伝子発現の比較には、要約統計量の計算、パターンの特定、データの分布の可視化がしばしば必要となります。Rパッケージはこれらのタスク専用に設計されており、基本的なコマンドラインツールの機能を超える機能を提供します。

7. 考慮事項とベストプラクティス

  • 適切な正規化手法の選択:

    • サンプル内比較: 同一サンプル内の異なる遺伝子間の相対的な発現量を比較するには、FPKMおよびTPMが使用できます 2。
    • サンプル間比較: サンプル間の比較には、TPMが一般的にFPKMよりも推奨されます。これは、TPMの合計がサンプル間で一定であるためです 2。CPMは、シーケンス深度で正規化されたカウントの基本的な比較に使用できますが、遺伝子長を考慮していません 1。
    • 差分発現解析: edgeRやDESeq2のようなツールでは、ローカウントを入力として使用することが一般的に推奨されます。これらのツールは、組成バイアスを考慮し、データの分布をモデル化するより高度な正規化手法（例えば、edgeRのTMM 1 やDESeq2のサイズ因子正規化 8）を採用しています。

  • 差分発現におけるFPKMとTPMの限界: これらの手法は主に正規化を目的としており、正規化の過程で分散に関する情報が失われる可能性があるため、差分発現の統計的検定には適していません 8。

  • 実験計画とバッチ効果補正の重要性: 適切な実験計画（十分なレプリケート数など）と潜在的なバッチ効果への対処は、RNA-Seq解析で信頼性の高い結果を得るために非常に重要です。これらの要因は、関心のある生物学的条件とは独立して遺伝子発現の測定に影響を与える可能性があります。

  • small RNA-Seqにおける考慮事項: small RNA-Seqでは、一部の高度に発現するsmall RNAが分布を歪める可能性があるため、TPMには限界があるかもしれません 17。

正規化手法の選択は、特定の研究課題と比較の目的に基づいて行うべきです。万能なアプローチは存在しません。異なる正規化手法は、異なる種類のバイアスに対処します。各手法の長所と短所を理解することは、特定の解析に最も適切な手法を選択するために不可欠です。例えば、サンプル間の全体的な転写活性を比較することが目的であれば、TPMが適しているかもしれません。個々の遺伝子の差分発現を特定することが目的であれば、カウントデータを直接モデル化する手法が推奨されます。

8. 結論

本レポートでは、RNA-Seqデータ解析における主要な正規化手法であるFPKM、TPM、およびCPMについて、その定義、計算方法、そして重要な差異を解説しました。FPKMは遺伝子長とシーケンス深度を考慮した正規化手法であり、主にペアエンドデータに用いられます。TPMも同様に両方の要素を考慮しますが、正規化の順序が異なり、サンプル間の比較においてより優れた性能を発揮します。CPMはシーケンス深度のみを正規化する最も基本的な手法です。

それぞれの正規化手法には適切な利用場面があり、特にTPMはサンプル間の比較においてその利点が強調されます。しかし、FPKMとTPMは差分発現解析には直接的には適しておらず、そのような解析にはローカウントと、より高度な統計モデルを用いた正規化手法が推奨されます。RNA-Seqデータの正確な解釈のためには、これらの正規化手法を理解し、実験の目的やデータの特性に応じて適切な手法を選択することが不可欠です。