- 著者: Jeff Jones, Elliot J. MacKrell, Ting-Yu Wang, Brett Lomenick, Michael L. Roukes, Tsui-Fen Chou
- Corresponding author: Jeff Jones (Proteome Exploration Laboratory, Beckman Institute, California Institute of Technology, Pasadena, CA)
- 雑誌: BMC bioinformatics
- 発行年: 2023
- Epub日: 2023-06-06
- Article種別: Original Article
- PMID: 37280522
背景
質量分析 (MS) ベースの定量プロテオミクスは、複雑な生物学的システムの解明やバイオマーカーの探索において極めて重要な役割を担う [1-7]。しかし、ProteomeDiscoverer、MaxQuant、Skyline、DIA-NN (data-independent acquisition neural networks) など、多様な測定・カタログ化プラットフォームが存在し、それぞれが独自のデータ処理、フィルタリング、変換、可視化手法を持つ [8-17]。この結果、出力フォーマットが非標準化され、データ構造や変数名が混在するため、研究者はデータ取り込み、クレンジング、解析のために場当たり的なスクリプトを開発せざるを得ず、保守性の低いコード環境を生み出す原因となっていた。例えば、Subramanian et al. ProcNatlAcadSciUSA 2005 や Kugeratski et al. NatCellBiol 2021 のような機能エンリッチメント解析は、標準化されたデータ入力形式を必要とするが、現状ではその準備に多大な労力を要する。
既存のRベースの後解析パッケージとして、pmartR [18]、protti [19]、DEqMS (differential expression analysis of MS data) [20] などが挙げられる。また、ウェブ実装を伴うMSstats (Mass Spectrometry Statistical Analysis) [21, 22]、DAPAR (differential analysis of protein abundance with R) [23, 24]、ProteoSign [25] といったツールも存在する。しかし、これらのツールは特定の解析手法セットに特化しているか、あるいは複数の解析プラットフォームからのデータ取り込みに十分に対応しておらず、プラットフォーム横断でエンドツーエンドに連結できる統一されたワークフロー基盤が不足していた [17, 26-35]。例えば、タンパク質カウントの棒グラフ作成はRのtidyverseパッケージで比較的容易だが、複数のサンプル間のタンパク質カウントの重複をVenn図で可視化するには、一貫性のない高度なデータ操作が必要となる。この点において、既存のツールではデータ探索の柔軟性や再現性の確保に課題が残されている。
さらに、正規化を補完の前に行うか後に行うか、あるいはヒトケラチンなどのコンタミネーションをいつ除去するかといった「ラストマイル解析」の選択肢を、研究者が自由に並べ替えて探索できる柔軟な環境も手薄であった。例えば、コンタミネーションがサンプル間で異なる場合、正規化前に除去することが推奨されるが [37]、共培養された生物のデータを除去する場合は正規化後が望ましい場合もある。このような解析上の考慮事項に対して、シンプルかつ容易な実装で探索できる環境が未確立であった。
本研究で開発されたtidyproteomicsは、このような知識のギャップを埋めることを目的とし、複数のプロテオミクスプラットフォームからのデータ探索、個別の関数と解析順序の制御、論理的なフローでの複雑なワークフロー構築を可能にするオープンソースRパッケージとして開発された。tidyverseの設計哲学（共通の文法、データ構造でデータ操作をエンドツーエンドに連結）に触発され、複雑な解析を小さなステップ単位に分解することで、データ解析の柔軟性と再現性を向上させることを目指した。
目的
本研究の目的は、ProteomeDiscoverer、MaxQuant、Skyline、DIA-NNといった複数の質量分析ベースの定量プロテオミクス解析プラットフォームからの出力を、単一の統一データオブジェクトに取り込み、データキュレーション、正規化、補完、差次的発現解析、機能エンリッチメント解析、および可視化を、tidyverseに準拠したパイプライン構文で実行できるオープンソースRパッケージ「tidyproteomics」を設計し、公開することである。
具体的には、以下の機能を提供することを目指した。
- 多様なプロテオミクスデータ出力形式を、シンプルで標準化されたデータオブジェクトに変換する機能。これにより、データ取り込みと前処理の複雑性を軽減する。
- データ操作、正規化、補完、差次的発現解析、機能エンリッチメント解析といった一連の解析ステップを、エンドツーエンドで連結可能なパイプラインとして提供する機能。これにより、ユーザーは複雑な解析ワークフローを簡潔に記述できる。
- 各解析関数において、ユーザーが任意の手法を選択したり、独自のアルゴリズムを組み込んだりできる柔軟性を提供し、解析選択が結果に与える影響を研究者自身が探索的に制御できる枠組みを構築すること。
- 解析の各ステップを自動的に追跡・記録し、文献引用情報とともに保存することで、解析の完全な再現性を確保するためのメカニズムを実装すること。
- プログラミング経験の少ないユーザーでも利用できるよう、ウェブベースのR Shinyアプリケーションを提供すること。これにより、ツールのアクセシビリティを向上させる。
これらの目的を達成することで、複雑なプロテオミクスデータ解析ワークフローを簡素化し、データ探索の柔軟性を高め、研究コミュニティにおけるプロテオミクスデータ解析の標準化と再現性向上に貢献することを目指す。
結果
統一データオブジェクトの4テーブル構造による効率的なデータ管理: import() 関数は、ProteomeDiscoverer、MaxQuant、Skyline、DIA-NNなどの多様な質量分析プラットフォームからの出力を、4つの非冗長テーブル (Experiments, Quantitative, Accounting, Annotations) に正規化する (Table 1)。Experimentsテーブルはsample_id、import_file、sample_file、sample、replicateといった実験メタデータを格納する。Quantitativeテーブルは、8文字のsample_id、identifier、およびraw値とmedian、linear、loess、randomforest正規化後のabundance値を保持する。Accountingテーブルは、補完率 (imputed) を0-1値で保持し、peptides、unique_peptides、proteinsの整数カウントを記録する。Annotationsテーブルは、GO (Gene Ontology)、KEGG (Kyoto Encyclopedia of Genes and Genomes)、UniProtなどの生物学的注釈を一対多構造で格納し、主要な定量データに干渉することなく追加できる。この「冗長情報を別テーブルに分離」する現代的データ構造哲学により、単一テーブルのサイズと複雑性が削減され、主要コンポーネントへのアクセス速度が向上する。collapse() 関数を使用することで、ペプチドデータをタンパク質レベルに集約でき、protein inferencingアルゴリズム、タンパク質数、要約関数 (sum、median、meanなど)、共有ペプチド按分の選択が可能である。
正規化の自動選択とrandomforestの優位性: HCT116 p97ノックダウン実証データセットを用いた解析では、タンパク質カウントがctrl群で5729 ± 59、kndw群で5817 ± 37 (各群n=3 biological replicates、95% CI表示、match-between-runsあり) と算出された。両群のタンパク質重複は、Venn図およびEuler図で6015共通、497群特異的として可視化された (Fig 2C, D)。normalize() 関数を用いて5種類の正規化手法を比較した結果、dynamic range (95%CI log10) とquantitative variationにおいて、rawデータがrange 5.6 / CV 0.21、median正規化がrange 5.7 / CV 0.25、limma正規化がrange 5.6 / CV 0.19、randomforest正規化がrange 5.5 / CV 0.17を示した (Fig 3B)。この結果から、randomforest正規化が全体のダイナミックレンジをほとんど損なわずに分散を最も大きく低下させることが示された。plot_dynamic_range() による密度ヒートマップ (Fig 3D) は、CVが低定量域で増大する傾向を示し、randomforest正規化がこの低存在量領域の大きな分散を最も良好に抑制することを実証した。select_normalization() 関数は、これらのCV、dynamic range、PCA成分変動を統合した加重スコアに基づいて、最適な正規化手法を自動推奨する。
差次的発現解析と新規可視化機能: expression() 関数は、2群間の比 (例: kndw/ctrl、各群n=3 replicates) として差次的発現を定義し、StudentのT検定またはlimmaパッケージのempirical Bayes法を選択できる。解析結果は、従来型のplot_volcano() (Fig 5A) および、定量的寄与度を強調する新規のplot_proportion() (Fig 5B) で可視化される。plot_volcano() はlog2_foldchange軸が約-1.6から+1.6、p=0.0001まで描画される。plot_proportion() は、プルダウン実験のように差が大多数に及ぶと予測される場面で特に有用であり、log2倍率1.6 (約3.0-fold) 級のup/down制御タンパク質 (RACGAP1 (Rac GTPase activating protein 1)、UBE2C (Ubiquitin conjugating enzyme E2 C)、SQLE (Squalene epoxidase)、DLG3 (Discs large homolog 3) など) を色分け表示する。さらに、2つの発現解析を対比する新機能plot_compexp() を導入し、同一controlに対する異なる処理 (別化合物・別遺伝子変異) や、Wilcoxon rank sum法とempirical Bayes法という異なる検定法の比較が可能になった (Additional file 1: Fig. S1)。enrichment() 関数は、GSEA [41] またはWilcoxon rank sum法で機能エンリッチメントを算出し、termをp value (0.25-1.00) で順位付けし、size 1000-3000のバブルとしてplot_enrichment() で描画する (Fig 5C)。
品質管理、解析追跡、およびアクセシビリティ:
summary() 関数は、subset() と同一のセマンティクスで任意の変数セットの要約統計を生成し、plot_counts()、plot_quantrank()、plot_venn()、plot_euler()、plot_normalization()、plot_variation()、plot_variation_pca()、plot_heatmap()、plot_pca() を含む包括的な可視化関数群を提供する (Fig 2, Fig 3, Fig 4)。Fig 4では、階層的クラスタリングヒートマップと第1-2主成分のPCAスコアプロットにより、ctrl群とkndw群のアンバイアスなクラスタリングが示された。subset() 関数は、%like% 正規表現演算子を用いて意味的フィルタリング (例: description %like% ‘ribosome’) を可能にし、merge() および reassign() 関数で複数ソースのデータ結合や群割り当てを行える。tidyproteomicsの最大の特徴の一つであるoperations() 関数は、各変換ステップを文献引用付きで自動記録し、解析の完全な再現性を担保する。プログラミング経験を問わないR Shinyウェブアプリの提供により、初学者から熟練Rユーザーまで幅広い層が利用できるようになった。
考察/結論
tidyproteomicsは、tidyverseの設計哲学（共通の文法・データ構造でデータ操作をエンドツーエンドに連結）をLCMS (liquid chromatography mass spectrometry; 液体クロマトグラフィー質量分析) プロテオミクス後解析へ適用したツールセットである。既報のRパッケージ (pmartR、protti、DEqMS、MSstatsなど) と異なり、本研究では単一の統一データオブジェクトを介して複数解析プラットフォーム (ProteomeDiscoverer、MaxQuant、Skyline、DIA-NN) の出力を横断的に扱い、関数の実行順序を研究者が自由に並べ替えられる点が対照的である。これにより、ユーザーは解析選択が結果に与える影響を探索的に制御できる。
plot_proportion() による定量的寄与度可視化、2つの解析結果を対比するplot_compexp()、および文献引用付きで解析ステップを自動追跡するoperations() 関数は、本パッケージで新規に導入された機能である。これらの機能は、これまで報告されていない再現性確保の枠組みを提供し、複雑なプロテオミクス解析の透明性と信頼性を向上させる。
臨床応用としては、エクソソームやマイクロベシクルなどの細胞外小胞 (EV) プロテオーム研究において、その有用性が期待される。例えば、複数患者検体からのLC-MS/MSデータを統一フォーマットで管理し、EV特異的バイオマーカーの差次的発現解析や機能エンリッチメント解析を、標準化された再現可能なワークフローで実施するためのbench-to-bedsideの橋渡しツールとして直接応用できる。これは、EV表面タンパク質の同定や、疾患診断・予後予測バイオマーカーの探索に貢献しうる。
本研究のlimitationとして、Tukey’s Median PolishやMaxLFQ [52] といった近年の高度なタンパク質定量アルゴリズム、およびPSM (peptide spectrum match; ペプチドスペクトルマッチ) カウントを考慮した検定 (DEqMS、MSqRob (mass spectrometry robust regression) 系) が未統合である点が挙げられる。これらのアルゴリズムの組み込みは、今後の課題となる。また、Table 2が提示する実験デザイン別の正規化/補完戦略 (例: 小群間変化なら任意正規化 + randomforest between-group補完、アフィニティ捕捉ならcommon contaminants除去 + minimum within-group補完など) は、著者の経験則に基づくものであり、比較検証を経ていない。したがって、ユーザー側には十分なドメイン知識が前提となる。
関連するEV定量プロテオミクス研究 (例: Choi et al. JExtracellVesicles 2012、Kugeratski et al. NatCellBiol 2021) や、同時期に登場したlabel-freeプロテオミクス向けRパッケージ (例: Ranathunge et al. BioinformAdv 2023) と並べて運用することで、EV表面プロテオーム (Doulabi et al. CommunBiol 2022) の標準化解析パイプラインを構築できる可能性を秘めている。
方法
本論文は新規アルゴリズムの提案ではなく、質量分析ベースの定量プロテオミクスデータ解析のためのソフトウェア基盤の記述を主とするSoftware Articleであるため、前向きコホート研究や動物実験は含まれない。RパッケージtidyproteomicsはMITライセンスのオープンソースとしてGitHub (https://github.com/jeffsocal/tidyproteomics) に公開されており、プラットフォーム非依存で動作し、その他の要件は特にない。
機能実証用データセットには、HCT116ヒト大腸がん細胞株を用い、wildtype (ctrl shRNA) とp97/VCP (valosin-containing protein) 遺伝子の単一標的ノックダウン (kndw, p97 shRNA) の生物学的反復 (各n=3 replicate) をProteomeDiscoverer 3.0で解析したlabel-freeデータを使用した [38]。このデータセットはパッケージ内にプリロードされており、全ての図を再現するためのRスクリプトがSupplemental Materialsとして提供されている。
データ構造とキュレーション:
import() 関数は、ProteomeDiscoverer、MaxQuant、Skyline、DIA-NNなどの多様なプラットフォームからの出力を、4つの非冗長テーブル (Experiments, Quantitative, Accounting, Annotations) に正規化する (Table 1)。これにより、データ構造の簡素化と主要コンポーネントへのアクセス速度向上が図られる。collapse() 関数は、ペプチドデータをタンパク質レベルに集約し、タンパク質推論アルゴリズム、タンパク質数、要約関数 (sum, median, meanなど)、共有ペプチド按分の選択を可能にする。subset() 関数は、tidyverseのdplyrパッケージのfilter関数と同様のセマンティックな式を用いてデータをフィルタリングできる。また、%like% 正規表現演算子を導入し、パターンマッチングによるデータサブセット化を可能にした。merge() および reassign() 関数は、複数ソースのデータ結合や群割り当てに用いられる。
正規化と補完: normalize() 関数は、median shift、linear、loess、limma、randomforestの5種類の正規化手法を一括実行できるラッパーとして設計された。select_normalization() 関数は、変動係数 (CV)、ダイナミックレンジ、主成分分析 (PCA) の上位3成分の変動を統合した加重スコアに基づいて、最適な正規化手法を自動選択する。impute() 関数は、欠損値補完をサポートし、missForestパッケージのrandom forest法 [47] および最小値補完を、サンプルグループ内またはグループ間で適用できる。これにより、MNAR (missing not at random; 非ランダム欠損) が想定される場合の挙動も検証可能である。
差次的発現解析と機能エンリッチメント: expression() 関数は、2群間の比 (例: kndw/ctrl) として差次的発現を定義し、StudentのT検定またはlimmaパッケージのempirical Bayes (経験ベイズ) 法を選択できる [54]。enrichment() 関数は、Subramanian et al. ProcNatlAcadSciUSA 2005 [41] またはWilcoxon rank sum (順位和) 検定を用いて機能エンリッチメントを算出する。
可視化と品質管理: plot_volcano()、plot_proportion()、plot_enrichment()、plot_heatmap()、plot_pca() など、多数のplot_() 関数群が可視化機能を提供する。plot_proportion() は、定量的寄与度を強調する新規の可視化手法である。plot_compexp() は、2つの発現解析結果を対比する新機能として導入された。summary() 関数は、subset() と同様のセマンティクスで任意の変数セットの要約統計を生成し、plot_counts()、plot_quantrank()、plot_venn()、plot_euler()、plot_normalization()、plot_variation()、plot_variation_pca() などの包括的な可視化関数群を提供する。
解析追跡とアクセシビリティ: operations() 関数は、各変換ステップを文献引用付きで自動記録し、解析の完全な再現性を担保する。プログラミング経験を問わないユーザー向けに、R Shinyウェブアプリ (http://bioinformatics.pel.caltech.edu/tidyproteomics/) も開発され、アクセシビリティが向上している。本研究では、統計手法としてStudentのT検定、empirical Bayes法、Wilcoxon rank sum検定が用いられている。