- 著者: Stefan Graw, Jillian Tang, Maroof K Zafar, Alicia K Byrd, Chris Bolden, Eric C. Peterson, Stephanie D Byrum
- Corresponding author: Stephanie D Byrum (Department of Biochemistry and Molecular Biology, University of Arkansas for Medical Sciences, Little Rock, AR; Arkansas Children’s Research Institute, Little Rock, AR)
- 雑誌: ACS Omega
- 発行年: 2020
- Epub日: 2020-09-30
- Article種別: Original Article
- PMID: 33073088
背景
質量分析 (mass spectrometry) 技術の急速な進歩により、1 回の実験で数千のタンパク質を複数サンプルにわたって高精度・高速に定量することが可能となった。TMT (tandem mass tag) マルチプレックスおよびラベルフリー定量法は、現代プロテオミクス研究の標準的アプローチとして広く普及している。しかし、プロテオミクスデータには試料調製・装置キャリブレーション・温度変化・未知源の非生物学的バイアスが混入しやすく、適切な正規化なしでは誤った結論に至るリスクが高い。これらのバイアスはバッチ効果とは異なり測定・制御が困難であり、統計解析時の事後補正が難しいという特性を持つ。
正規化はこのようなシステマティックバイアスを補正しサンプル間の比較可能性を向上させる重要な前処理ステップであり、その多くは DNA マイクロアレイ技術から採用されてきた。Välikangas et al. (2018) による系統的評価では、VSN (variance stabilizing normalization)・線形回帰正規化・局所回帰正規化が一貫してトップランクの性能を示す一方、各データセットに最適な正規化手法は異なるという重要な知見が示された。正規化ツールとして、Chawade et al. (2014) は Normalyzer を開発し、PCV (pooled coefficient of variation)・PMAD (pooled median absolute deviation)・PEV (pooled estimate of variance) による複数正規化手法の評価機能を提供した。Webb-Robertson et al. (2011) は 8 種のペプチド選択法と 5 種の正規化法を組み合わせた評価ツール SPAN を提供した。
しかし、これらの先行ツールには共通して大きな課題があった。正規化評価から欠損値補完 (missing value imputation)・差次発現解析 (differential abundance analysis) 手法の比較・統計パワー計算まで、プロテオミクス解析の全ワークフローを一貫して提供する公開ツールが不足していた。さらに R 言語に不慣れな研究者が GUI (graphical user interface) だけで全ステップを完結できる環境も手薄であり、解析手法の選択と標準化における知識のギャップが明確に存在していた。この gap in knowledge を埋め、研究者がデータセット固有の最適正規化手法を容易に評価・選択できるツールの開発が強く求められていた。
目的
TMT6plex・TMT10plex マルチプレックスおよびラベルフリー質量分析プロテオミクスデータに対して、正規化手法の体系的評価・欠損値補完・差次発現解析手法の比較・統計パワー計算を一貫した GUI ワークフローで実行できる R-Shiny ツール「proteiNorm」を開発・公開すること。R 言語に不慣れな研究者でも MaxQuant 出力ファイルから点とクリックの操作で全解析ステップを実行できる環境を提供する。加えて、弱シグナル・強シグナル・既知スパイクインという性質の異なる 3 種の代表的プロテオミクスデータセットを用いて proteiNorm の有効性を実証し、実験デザインに応じた最適正規化戦略の選択基準を明確にすることを目的とする。
結果
proteiNorm のアーキテクチャと主要機能: proteiNorm は R Shiny を基盤とした GUI ツールとして実装され、ポイントアンドクリックで全解析ステップを実行できる (Figure 1)。ワークフローは (1) MaxQuant 出力の peptides.txt 入力・Top3 フィルタリング (任意)、(2) proteinGroups.txt 入力・メタデータ設定 (サンプル名・グループ・バッチラベル)、(3) 外れ値・不要サンプル除外 (「filter」タブ)、(4) 8 種の正規化手法の性能比較と適切手法の選択 (「Normalization」タブ)、(5) 欠損値パターン可視化と補完手法の選択 (任意)、(6) 「DAtest」タブで差次発現解析手法の自動比較・ランク付けと正規化済みデータのエクスポート、(7) 「DApower」タブで効果量範囲にわたる統計パワー推定、という 7 ステップで構成される。各正規化手法の性能は総強度分布 (Figure 2A)・PCA プロット (Figure 2B)・欠損値ヒートマップ (Figure 2C)・PCV/PEV/PMAD (Figure 2D-F)・グループ内相関 (Figure 2G)・相関ヒートマップ (Figure 2H)・log2 比率密度分布 (Figure 2I) の 9 種の可視化で自動比較される。ツールは「Save Figures」ボタンで全図の保存が可能であり、GitHub と Shiny Web アプリの両形態で公開されている。
マウス TMT6plex データ (弱シグナル) における正規化評価: マウス脳線条体のメタンフェタミン処置 vs 食塩水対照の TMT6plex データを解析した。メタデータはバッチ 1-2、nucleus accumbens (NA) と striatum (S)、処置 (T) と対照グループに分類されている (Table 1)。解析は striatum サンプル n=6 (処置群 n=3、対照群 n=3) に限定し、nucleus accumbens サンプル n=6 は「filter」タブで除外した。いずれのサンプルも総強度分布・PCA に基づく外れ値は認められなかった (Figure S1-S2)。正規化手法の比較では、VSN がグループ内変動の 3 指標すべてで最良の性能を示した (Figure 2D-F)。PCV では VSN が約 0.012 を達成したのに対し、median・mean・quantile 等の他手法は 0.013-0.014 の範囲にとどまった。log2 比率密度分布では VSN と RLR のみがゼロ周辺に中心化され、他の正規化手法は程度の差こそあれ偏りを残した (Figure 2I)。欠損値ヒートマップでは明確な MAR/MNAR パターンは認められなかった (Figure 2C)。limma (Ritchie et al. 2015; Ritchie et al. NucleicAcidsRes 2015) は完全データを必要としないため imputation なしで解析を継続した。PCA では薬物処置効果よりも試料個体差 (nucleus accumbens vs striatum) が主成分の分散を支配しており、striatum の処置・対照グループ間の明確な分離は確認されなかった (Figure 2B)。弱シグナルデータでは正規化手法の選択が最終的な有意タンパク質リストに影響を与える可能性が示唆された。
乳癌細胞株 TMT10plex データ (強シグナル) における正規化評価: MCF10A (非腫瘍性上皮細胞株)・MCF7 (ER/PR+ 乳癌)・HCC1954 (HER2+ 乳癌) の 3 細胞株に対して hydroxyurea 処置・無処置の比較を TMT10plex で実施した (n=18 サンプル + pooled n=2、バッチ 3-4、Table 2)。Pooled サンプルはバッチ効果評価後に解析から除外した。欠損値ヒートマップでは、欠損がバッチおよび細胞株-処置の組み合わせと相関する MNAR パターンが観察された (Figure S8C)。バッチ効果は raw peptide データで顕著であり、filtered protein データでは縮小することが PCA により確認された (Figure S7)。cyclic loess 正規化を選択し、全正規化手法で総強度分布はほぼ均一化された (Figure S8A)。PCV・PEV・PMAD の 3 指標では cyclic loess が他の手法と同等またはわずかに優れた性能を示した (PCV: cyclic loess 約 0.013 vs median 約 0.014) (Figure S8D-F)。PCA では正規化後に細胞株ごとのクラスター分離が明確となった (Figure S8B)。相関ヒートマップでは、同一細胞株の処置・無処置サンプル間の Pearson 相関が最も高く (r > 0.98)、異なる細胞株・異なる処置サンプル間 (異バッチ) の相関が最低値を示した (Figure S8H)。強シグナルデータでは複数の正規化手法が同等の性能を示し、手法選択の影響が弱シグナルデータと比べて相対的に小さかった。
スパイクインデータ (ラベルフリー) における正規化性能の客観的検証: Yeast lysate に UPS1 標準タンパク質を既知濃度 (50, 125, 250, 500, 2,500, 5,000, 12,500, 25,000, 50,000 amol) でスパイクした公開データ (ProteomeXchange ID: PXD001819、各濃度 n=3 replicate、Table 3) を解析した。濃度 50-500 amol では測定値がほとんど得られなかったため解析から除外した。各濃度グループの replicate 平均を用いて異なる濃度間の比率を算出し、raw データと VSN 正規化後データを比較した (Figure 3)。VSN 正規化後データは raw データより一貫して真の濃度比率 (灰色破線) に近い推定値をもたらした。2,500/50,000 amol の比率 (理論値 0.05) では、raw データの推定値が約 0.08 であったのに対し、VSN 正規化後は約 0.04 に改善した。25,000/50,000 amol の比率でも正規化後の精度向上が明確であった。このスパイクイン実験により、スパイクイン標準を用いた客観的評価基準のもとでも VSN 正規化が真のシグナルをより正確に推定できることが定量的に実証された。特に、スパイクイン濃度と強度が高い領域ほど正規化の効果が顕著であることが明らかとなった。
差次発現解析手法の比較と統計パワー評価 (DAtest 統合): 正規化後データに対して DAtest パッケージ経由で複数の統計検定手法の性能を自動比較する機能を実装した。DAtest はデータをシャッフルし指定した効果量 (例: 2-fold change に相当する spike-in parameter = 2) でスパイクして、各手法がスパイクインした特徴を正しく同定できるか・FDR を適切に制御できるかを繰り返し評価する。データタイプ (カウントデータか強度データか) に応じて適切な統計手法が自動的に選択候補として提示される。結果はランキング表と図で出力され、研究者は最適な差次発現解析手法を客観的データに基づいて選択できる。DApower タブでは、指定した fold change 近傍の効果量範囲 (例: 1.5-3x-fold) にわたる統計パワー (真のシグナルを正しく検出する確率) を各手法で推定し、研究デザインの適切性を事前評価できる。計算コストの高い permutation ベースの統計手法は除外オプションにより実行時間を大幅短縮できる。正規化済みデータは「DAtest」タブの「Save file」ボタンからエクスポートして、limma 等の下流ツールでの差次発現解析に使用できる。
考察/結論
先行ツールとの違い: SPAN (Webb-Robertson et al. 2011) や Normalyzer (Chawade et al. 2014) とは異なり、proteiNorm は正規化評価に留まらず、欠損値補完手法の可視化と選択・複数の差次発現解析手法の自動比較・統計パワー計算までの一連のプロテオミクス解析ワークフローを一つの GUI で提供する。これまでの研究では、正規化ステップと下流の差次発現解析を橋渡しする統合的なツールは対照的に存在しておらず、研究者が複数ツール間でデータを受け渡す必要があった。また、RLR の実装は Chawade らの Normalyzer に倣っているが、既報の Normalyzer は TMT マルチプレックスデータと DAtest による差次発現解析比較機能を統合していない。これらの相違が proteiNorm の主要な差別化要因である。
新規性: 本研究で初めて、弱シグナル (マウス脳)・強シグナル (乳癌細胞株)・既知スパイクインという性質の異なる 3 種の実験デザインを代表するプロテオミクスデータセットを同一ツールで系統的に評価し、最適正規化手法が実験デザインに依存することを定量的に実証した。マウスデータ (弱シグナル) では VSN が全 3 グループ内変動指標で最良を達成し、乳癌データ (強シグナル) では cyclic loess が最適として選択された。スパイクインデータでは VSN による正規化が真の濃度比率推定の精度を改善した (2,500/50,000 amol で raw 0.08 → 正規化後 0.04)。「VSN と cyclic loess は一貫してトップランク」という知見はVälikangas et al. (2018) の先行研究と整合しており、本研究で初めてユーザーフレンドリーな GUI ツールとして実装・公開したことが新規の貢献である。また、DAtest を統合することで差次発現解析手法の体系的比較を同一プラットフォーム内で実施できる機能は、これまで報告されていない統合的アプローチである。
臨床応用と研究への意義: 不適切な正規化による誤った結論は、バイオマーカー発見や創薬ターゲット探索において多大なコスト・時間の浪費につながる。proteiNorm の開発はプロテオミクス研究の再現性・標準化向上に直接的な臨床的意義を持つ。特に、腫瘍プロテオミクス研究 (エクソソーム・細胞外小胞のプロテオーム解析含む) では TMT によるマルチサンプル・マルチバッチ比較が標準的に行われており、本ツールはそのような解析の前処理標準化に具体的に貢献できる。R 言語スキルを必要としない GUI 設計は、臨床医・生物学者が直接データ解析に参加できる体制を構築するための bench-to-bedside のブリッジとして機能し、学際的プロテオミクス研究の推進に資する。Välikangas et al. による系統的評価が示すように、大多数の正規化手法は単純な log2 変換より優れており、本ツールが最適手法選択の敷居を下げることで研究全体の信頼性が向上する。
残された課題と今後の検討: 現バージョンの proteiNorm にはいくつかの limitation がある。第 1 に、スパイクイン標準やプロテオミクスリファレンスを用いた正規化評価が現実装では未統合であり、スパイクインを設計していない通常実験では客観的な評価基準を用いた比較が困難なままである。第 2 に、現在の入力フォーマットが MaxQuant 出力に特化しており、DIA (data-independent acquisition) を解析する DIA-NN や Spectronaut、あるいは Proteome Discoverer 等の他のプラットフォームへの対応拡張が今後の課題として残されている。第 3 に、バッチ効果補正 (ComBat 等) が正規化手法として含まれていない点も limitation であり、マルチバッチデザインでは別途バッチ補正が必要な場合がある。将来的には、スパイクイン評価機能の統合・追加の正規化・imputation 手法の実装・入力フォーマットの拡張が予定されており、更なる検討を要する。大規模コホート (n=100 以上) での計算効率や単一細胞プロテオミクスへの適用可能性についても今後の研究が必要である。
方法
ツールの実装と公開: proteiNorm は R 言語 (version 3.6) と Shiny パッケージを用いて実装され、Web アプリケーション (https://sbyrum.shinyapps.io/proteiNorm/) および GitHub リポジトリ (ByrumLab/proteiNorm) として無償公開されている。ワークフローは「data」→「filter」→「Normalization」→「DAtest」→「DApower」の 5 タブ構成で、データ入力から正規化評価・欠損値補完・差次発現解析・統計パワー評価まで順に進む設計となっている。
データ入力形式: MaxQuant 出力の peptides.txt (任意) と proteinGroups.txt を直接入力でき、TMT 実験では「Reporter intensity corrected」で始まるカラム名、ラベルフリー実験では「Intensity」で始まるカラム名を自動認識する。ツールはデータタイプを自動判定し、MaxQuant が common contaminants および reverse sequences にフラグを立てたタンパク質を自動除外する。TMT 実験では Thermo Fisher のロット固有の同位体分布補正係数 (-2, -1, +1, +2 carbon isotopes) が適用された「Reporter intensity corrected」値を使用する。
ペプチドフィルタリング: Top3 法 (Silva et al. 2006) により、各タンパク質の最も強度の高い 3 ペプチドの強度を平均してタンパク質強度を算出し、外れ値ペプチドによる相対存在量の歪みを低減したフィルタリング済み proteinGroups ファイルを生成する (任意ステップ)。
実装した正規化手法: log2 変換、median 正規化、mean 正規化、VSN (Huber et al. 2002)、quantile 正規化 (Bolstad 2019)、cyclic loess 正規化、RLR (global robust linear regression, Chawade et al. 2014)、global intensity 正規化 の計 8 種類を実装した。これらは DNA マイクロアレイおよびプロテオミクスで広く用いられてきた手法群である。
性能評価指標: PCV (pooled intragroup coefficient of variation)、PEV (pooled intragroup estimate of variance)、PMAD (pooled intragroup median absolute deviation)、グループ内 Pearson 相関係数、log2 比率密度分布 (全 2 群組み合わせ) の 5 種類の定量的指標と、総強度分布・PCA (principal component analysis)・相関ヒートマップの可視化により各手法を多角的に評価する。グループ内変動指標は値が小さいほど良好で、相関指標は値が大きいほど良好な正規化を示す。
欠損値補完手法: KNN (k-nearest neighbors, Hastie et al. 2018)、QRILC (Quantile Regression Imputation of Left-Censored, Lazar 2015)、MinDet (deterministic minimal value imputation)、MinProb (stochastic minimal value imputation)、Min (minimal value imputation)、Zero の 6 種類を実装した。欠損値パターンの MAR (missing at random) / MNAR (missing not at random) の判別には DEP Bioconductor パッケージ由来の修正ヒートマップを使用する。
差次発現解析の比較: R パッケージ DAtest (Russel et al. 2018) を統合し、edgeR・limma 等に相当する複数の統計手法の FDR (false discovery rate) 制御能とスパイクイン特徴の検出感度を自動比較・ランク付けする。DAtest では完全データ (欠損なし) が必要であるため、適切な imputation 手法を選択した上で差次発現解析手法の評価を行う。DApower タブでは指定した fold change 近傍の効果量範囲にわたる統計パワーを各手法で推定できる。
検証データセット: (1) マウス線条体 TMT6plex データ (メタンフェタミン処置 vs 食塩水対照、striatum サンプル n=6、batchは 2 バッチ、弱シグナル)、(2) 乳癌細胞株 TMT10plex データ (MCF10A 非腫瘍性上皮細胞株・MCF7 ER/PR+ 細胞・HCC1954 HER2+ 細胞、hydroxyurea 処置・無処置、n=18 サンプル + pooled サンプル n=2、バッチ 3-4、強シグナル)、(3) 公開スパイクインデータ (yeast lysate + UPS1 標準タンパク質 50-50,000 amol スパイク、n=27 サンプル、ラベルフリー、ProteomeXchange ID: PXD001819) の 3 種を使用した。全 mass spectrometry データは ProteomeXchange ID: PXD018152 で公開されている。