RNA-seqの落とし穴
Pitfalls in RNA-seq
東京農業大・農学部
植物育種学研究室 石森 元幸
日本育種学会第148回講演会ワークショップ
データ解析の「落とし穴」
〜あなたのその解析、大丈夫ですか?〜
2025年9月10日(水)
札幌コンベンションセンター
(北海道大学大会)
https://www.nodai.ac.jp/campus/map/
自己紹介
Name: 石森 元幸 Motoyuki Ishimori
Age: 41
Specialty: 植物育種学・園芸学
バイオインフォマティクス等
Target Plants:
主に野菜(ナス・トウガラシ)
花き(トルコギキョウ・ハナスベリヒユ等)
研究経歴
2004-2008 農学部(東京農工大学・植物遺伝育種学研究室・平田豊教授)
2008-2010 修士課程(東京農工大学大学院・同上)
2011-2015 博士課程(東京大学大学院・園芸学研究室)
2015-2021 特任研究員(東京大学大学院・生物測定学研究室)
2021-2025 助教(東京大学大学院・園芸学研究室)
2025- 准教授(東京農業大学・植物育種学研究室)
DEG検出のためのRNA-seqの手順
1.実験全体のデザイン(一番重要)
2.シーケンシング条件
3.クオリティチェック(QC)
4.リードマッピング(リファレンスありの場合)
DEG: Differentially Expressed Gene(発現変動遺伝子)
5.カウントデータ取得
6.データ正規化
7.DEG検出(検定)
RNA-seqにおけるデータ正規化とは?
*Rパッケージ「TCC」のサンプルデータを使用
| G1_rep1 | G1_rep2 | G1_rep3 | G2_rep1 | G2_rep2 | G2_rep3 |
gene_1 | 34 | 45 | 122 | 16 | 14 | 29 |
gene_2 | 358 | 388 | 22 | 36 | 25 | 68 |
gene_3 | 1144 | 919 | 990 | 374 | 480 | 239 |
gene_4 | 0 | 0 | 44 | 18 | 0 | 0 |
gene_5 | 98 | 48 | 17 | 1 | 8 | 5 |
gene_6 | 296 | 282 | 216 | 86 | 62 | 69 |
RNA-seqのカウントデータ(行:遺伝子、列:サンプル)
遺伝子3、および総カウント数の違い
DEG検出の大前提:『適切』な正規化
二群間比較におけるのDEG検出(検定)の原理
(帰無仮説)
グループ1の平均 = グループ2の平均
(対立仮説)
グループ1の平均 ≠ グループ2の平均
!! DEG !!
非DEG
qRT-PCR(相対定量)における正規化(補正)
RNA-seqにおける正規化
サンプル
A
サンプル
B
増幅量(サイクル数)
調べたい遺伝子
サンプル
A
サンプル
B
増幅量(サイクル数)
ハウスキーピング
遺伝子
2倍
サンプル
A
サンプル
B
相対発現量
1/2倍
調べたい遺伝子
非DEG≒ハウスキーピング遺伝子を探す
| 正規化係数 |
G1_rep1 | 0.880 |
G1_rep2 | 0.860 |
G1_rep3 | 0.843 |
G2_rep1 | 1.085 |
G2_rep2 | 1.144 |
G2_rep3 | 1.188 |
TMM
(Trimmed Mean of M-values)
edgeRで採用
Robinson et al., 2010
主要な正規化手法
Size Factor Normalization
(median-of-ratios)
DESeq2で採用
Anders and Huber, 2010; Love et al., 2014
最も典型的な発現パターン*の基準サンプルを選出
基準サンプルに対する
M値(対数比)とA値(平均対数)を計算
極端なM値・A値の遺伝子*を除外
残った遺伝子のM値の加重平均から
正規化係数を算出
各遺伝子について全サンプルの幾何平均を計算
各サンプルについて各遺伝子のカウント数を
上記の幾何平均で割る
各サンプルの上記の幾何平均比の中央値を
Size Factorとする
Size Factorから
正規化係数を算出
・多くの場合で極端な違いは生じない
・高発現なDEGがある場合はDESeq2正規化が頑健
・多数のDEGがある場合はTMMがやや頑健
・DEGが多数ある場合は両者を比較すると良い
*TMMではこれらのパラメータを変更できることに注意
TMM/DESeq2正規化の比較
サンプル | 正規化係数(TMM) | 正規化係数(DESeq2) |
leaf_1 | 1.118 | 1.107 |
leaf_2 | 1.030 | 1.014 |
leaf_3 | 0.952 | 0.943 |
panicle_1 | 0.939 | 0.954 |
panicle_2 | 0.995 | 1.012 |
panicle_3 | 0.986 | 0.987 |
正規化手法の比較
(Li et al., PLOS ONE, 2017)
・DEGの存在自体が正しい正規化を阻む
(Kadota et al., 2012; 門田 2014)
・DEG同定 ⇒ 残りの非DEGのみで正規化
を繰り返すiDEGES(TCCパッケージ)も有効
(Sun et al., 2013)
← ソルガムRNA-seqデータ(DRA004664)
における正規化の例
edgeRとDESeq2のDEG検定
負の二項分布(negative binomial distribution)を仮定
Xij ~ NB(μij, φi)
Xij : サンプルjにおける遺伝子iのカウント, μij: 期待値(正規化後の発現量), φi: 遺伝子iの過分散パラメータ
RNA-seqのカウントデータ(生物学的反復)は
ポアソン分布では説明できない
過分散(over-dispersion)の性質を有する
一般化線形モデル(GLM)
線形予測子に群や処理条件などを組み込んだモデル
Dispersionの推定
DESeq2:遺伝子ごとに推定⇒全体トレンドに基づいて収縮
edgeR:共通・トレンド⇒経験ベイズで各遺伝子の分散推定
検定
DESeq2:Wald検定(発現変化がゼロか否か)、LRT
edgeR:準尤度(QL)F検定、LRT、exact test
両者の乖離が大きくなる条件
・反復数が最低限(3以下) ・Dispersion推定に差がある
・非常に高(低)発現遺伝子が多い
・デフォルトの閾値条件(edgeRはfold changeを考慮せず)
edgeRとDESeq2のDEG検出の比較
edgeRとDESeq2の比較1
(ソルガム公共データ、TCCパッケージ内の初期条件で実行)
edgeRとDESeq2の比較2
(正規化はTMM、検定は各初期条件で実行)
条件1でedgeRのみでDEG
条件1でDESeq2のみでDEG
非DEG用のカウントデータのスケーリング
RPKM
(Reads Per Kilobase per Million
mapped reads)
FPKM
(Fragment Per Kilobase per Million
mapped reads)
TPM
(Transcripts Per Million)
現在の標準的なスケーリング手法
遺伝子長で補正しないCPMが一般的
RPK(CPK)
(Reads/Counts Per Kilobase)
配列長のみの補正
RPK = カウント数×1000/配列長
RPKM = カウント数×1000/配列長×1000000/総リード数
FPKM = ペアエンド数×1000/配列長×1000000/総ペアエンド数
TPM =RPKM/総RPKM
・代表的なRNA-seqカウントデータのスケーリング方法
カウントデータのヒートマップ
(VST変換 ⇒ Zスコアスケーリング)
RPKM/TPM等がDEG検出向きでない理由
TMM (edgeR)とTPMの比較
(検定はedgeRのexact test)
TPM(上)のみでDEGの例
(検定はedgeRのexact test)
多群(3グループ以上),もしくは多因子(時間・処理・系統などの因子が3以上)から構成される多検体サンプルにおいては,以下の要因から正規化が難しくなる.
多群や多因子間の比較における正規化
多群間比較を必要とする実験デザイン
1.処理・時系列・系統などが複数ある多因子実験
2.一因子だが多群(例えば、時系列サンプルや処理区が3以上ある)
3.交互作用(例えば、処理×系統)を解析したい場合
*それぞれで二群間比較を行うとFDRの制御が難しくなり、誤った結論になることも
← 超多検体ではTMM正規化は不安定
(Lin et al., 2016)
きちんと正規化できているかを
確認することが重要
(基本的にはDESeq2が推奨)
多群や多因子比較におけるDEG検出
多群間比較におけるのDEG検出の基本原理
(帰無仮説)
全ての群で発現平均が等しい
(対立仮説)
少なくとも1つの群で発現平均が異なる
!! DEG !!
非DEG
二群および多群の場合の正規化
二群および多群におけるDEG検出
二群ではDEG、多群では非DEG
二群では非DEG、多群ではDEG
対応策:複数手法を出来るだけ比較する
DEG検出における4手法の比較
*各手法(パッケージ)のデフォルトの設定で解析した結果
・反復が多い場合はVoomかSAMseq (Soneson and Delorenzi, 2013)
・反復が少ない時はVoomかDESeq (Seyednasrollah et al., 2013)
・基本的にedgeRとDESeq2が最良(Ching et al., 2014)
・多群間比較ではedgeRかDESeq2 (Tang et al., 2015)
・NOISeqが最良で次点がedgeR (Stupnikov et al., 2021)
・EBSeqが最良だが反復が多ければDESeq2 (Li et al., 2022)
実験デザイン(反復数等)や目的(検出力重視か)によって「ベター」な手法は異なる
パッケージ | 分布モデル | 正規化方法 | 検定方法 | 特徴 | 適用例 |
edgeR | 負の二項分布 | TMM | Exact test/GLM-LRT GLM-QLF | 分散の経験ベイズ推定 ロバストオプション 複雑なデザイン、高速 | 小規模〜大規模 複雑デザイン |
DESeq2 | 負の二項分布 | median-of-ratios | Wald/LRT | 分散収束・収縮推定 log2FC収縮、再現性・頑健性 | 標準的な研究全般再現性重視 |
limma-voom | 正規線形モデル | TMM等 | Empirical Bayes (moderated t/F) | 平均–分散関係を精度重みに 変換高速・柔軟、大規模 | 大規模・反復の多い実験、多群デザイン |
baySeq | 負の二項分布 | 内部推定/オフセット | 経験ベイズ サンプリング | 群構成を柔軟に定義 計算負荷高 | 中〜大規模 複数群パターンの比較 |
EBSeq | 負の二項分布 | サイズファクター | Empirical Bayes | 状態確率で評価 | アイソフォーム解析 状態推定を重視 |
SAMseq | ノンパラメトリック | ライブラリサイズ補正置換再標本化 | ランク統計+置換(SAM) | 外れ値に強い 小標本では検出力低、計算量大 | 分布仮定なし |
NOISeq | ノンパラメトリック | CPM / RPKM / TMM | データ駆動型 ノイズモデル | 再現性指標 | 探索的解析・リード数不均一・バイアス |
まとめ:RNA-seq(DEG検出)の留意点
1.実験全体のデザイン(一番重要)
2.シーケンシング条件
3.クオリティチェック(QC)
4.リードマッピング(リファレンスありの場合)
5.カウントデータ取得
6.データ正規化
7.DEG検出(検定)
4*.トランスクリプトームアセンブリ
1*.ロングリード or マイクロアレイ
2*.シングル or ペアエンド・ストランド情報など
X.エンリッチメント解析(GO解析)
Y.次元圧縮と可視化(t-SNE or UMAP)
Z.クラスタリング
A.新世界~1細胞・空間トランスクリプトーム等
・実験デザイン~サンプリング~ライブラリ調製からデータ解析まで選択肢が多い
⇒ バッチ間の差異や手順・条件の違いによるアーティファクトへの注意が重要
・RNA-seqのみで重要な結論を得られることは意外なほど少ない(信頼性は高くない)
⇒ 研究の第一段階であり、qPCRやゲノム編集・組換え体によるvalidationが必要
・RNA-seqの有用性向上には非モデル植物における機能的アノテーションの充実が必要
⇒ 実はこれまで以上にゲノム編集などを活用した機能解析が重要になってくる?