#1．書こうと思ったいきさつ
以前こんなことがありました。
このツイートを見たとき、化合物の機序…は分からずとも、どのパスウェイにいればいいのかがわかれば部分的に達成できるかなあと思いMSEAを紹介してみました。ケモインフォのみなさま。gene ontology などを使った Gene set enrichment 解析のように、化合物セットを入力に、化合物の機序などの特徴をあぶりだす解析ツールってありますか? gene ontology のように化合物の定番の ontology ってありますか?
— ゲノムのほうの愛ちゃん🌙 🐧 (@dritoshi) November 1, 2019
普段僕は代謝物データの解析も行っているためこのツールには割と慣れていたんですが、意外とこの手の代謝物データ解析に関するツールってあまり知られていないのかなあと思いまして。化合物セットがどの代謝パスウェイに多く含まれているかを探すGSEAの代謝物版のMSEAっていうのがあるんですが、ちょっと想定されているものとは違いますかね…？ https://t.co/K5D3byhl39
— カミケ (@KKami1115) November 1, 2019
代謝物解析に関するお話ってとっつきにくい面があるのか…？と思ったので、今回はMetabolite Set Enrichment Analysis (MSEA)とそれを行うツール (サービス？)であるMetaboAnalyst、RパッケージであるMetaboAnalystRについての紹介です。
#2．MSEAについて
https://www.ncbi.nlm.nih.gov/pubmed/20457745
初出の論文はこれです。
GSEAと同様に、Metabolite Setと呼ばれる代謝物状態のライブラリを準備しておき、代謝物名リストや代謝物群の濃度をもとにエンリッチメント解析を行います。
使用するMetabolite Setですが、Pathway-basedなもの、Disease-basedなもの、Drug-basedなものと用途ごとに分かれています。最初のツイートだとパスウェイに限定してしまっていますが、パスウェイ以外の情報を使うことも可能です。
目的に応じて使い分けましょう。
アブストにもある通り、MSEAにはORA、SSP、QEAの3つの解析手法があります。前述の@dritoshi先生のように代謝物の名称のみが与えられる場合はORA (Overrepresentation Analysis)ですが、今回はQEA (Quantitative Enrichment Analysis)を行ってみましょう。
量的データとPathway-basedなMetabolite Setを使って、二群間比較において存在量が変動した代謝物群がどのパスウェイに多く含まれるのかを明らかにします。
#3．MetaboAnalystについて
https://www.metaboanalyst.ca/
上記手法を提唱したXia Labにて管理されている、種々の代謝物解析を行うためのWebサービスです。
今回行うMSEA以外にも様々な解析手法に対応しています。遺伝子・代謝物の統合解析1なんてのもできます。
#4．使用するデータとデータ加工
流石に生データからピーク表にもっていく方法を説明していると長くなってしまうので、色々流用します。
今回使用するデータはFukushima, et al. (2011)2にて取得された、Arabidopsis thaliana (シロイヌナズナ)のCol-0 (野生型)とtt4変異体、mto1変異体におけるGC-MSによる質量分析の結果です。
mto1変異体はCol-0と比較してメチオニンの過剰蓄積が起こるらしいので、MSEAが適切にそれを追えるか試してみましょう。
MetaboLightsという代謝物データ用DBの該当ページ3のFilesから、m_gc0002_metabolite profiling_mass spectrometry_v2_maf.tsvとs_GC0002.txtをダウンロード。
以下、MetaboAnalystに流せる形にデータを加工します。
library(tidyverse)
#データ読み込み
GC_table_raw = read_tsv("m_gc0002_metabolite profiling_mass spectrometry_v2_maf.tsv")
sample = read_delim("s_GC0002.txt", delim = "\t")
#アノテーションがついたピークのみに限定
GC_table = GC_table_raw[GC_table_raw$database_identifier!="unknown",]
#"Glycine (2TMS)"を削除し"Glycine (3TMS)"のみに (一般的に誘導体化は数が大きい方が優先されます)
GC_table = GC_table[GC_table$metabolite_identification != "Glycine (2TMS)",]
代謝物名について、そのままだとMetaboAnalystが認識しないためいろいろと調整 (バッドノウハウともいう)
GC_table$metabolite_identification = GC_table$metabolite_identification %>%
str_remove(pattern ="[:punct:][:digit:]*TMS[:punct:]") %>%
str_remove(pattern ="[:punct:].?MeOX[:punct:]") %>%
str_remove(pattern ="DL-") %>%
str_remove(pattern =fixed(" (1MEOX) ")) %>%
str_remove(pattern =fixed("D-(+)-")) %>%
str_remove(pattern =fixed("D-(-)-")) %>%
str_remove(pattern =fixed(" (1MEOX) (4TMS)")) %>%
str_remove(pattern =fixed("-2,2,3,3-d4")) %>%
str_remove(pattern =fixed("-3,4,5,6-d4")) %>%
str_remove(pattern =fixed("-13C5,15N")) %>%
str_trim(side = "right")
※MetaboAnalystに代謝物の名称を渡す際ですが、
| 表記 | 可否 |
|---|---|
| D-Glucose | Ok |
| D(-)Fructose | アカン |
| D(+)-Mannose | Ok |
| D-(+)-Glucose | アカン |
| Glucose | D-GlucoseとしてOk |
| このように、命名規則がバラバラでもMetaboAnalystは割と解釈してくれますが限度があります。 | |
| 後述するMetaboAnalystでの処理である程度カバーできますが、命名規則が統一されているかどうか確認しておきましょう。 |
Col-0とmto1の二群を残して、MetaboAnalystが指定する形4に整形します。
GC_table_me = GC_table[,22:75]
colnames(GC_table_me) = paste(sample$`Factor Value[Genotype]`, c(c(1:17), c(1:20), c(1:16)), sep = "_")
GC_table_me = GC_table_me[,c(str_subset(colnames(GC_table_me),"mto1"), str_subset(colnames(GC_table_me),"Col-0"))]
GC_table_me = rbind(as.numeric(as.factor(sample$`Factor Value[Genotype]`)), GC_table_me)
rownames(GC_table_me) = c("Label", GC_table$metabolite_identification)
GC_table_me = rownames_to_column(GC_table_me)
write_csv(GC_table_me, path = "gctableme.csv")
83種の代謝物と、Col-0が17サンプル、mto1が16サンプルある行列の完成です。
#5．MetaboAnalystの使い方
生成したgctableme.csvをhttps://www.metaboanalyst.ca/MetaboAnalyst/ModuleView.xhtml のEnrichment Analysisにアップロード。あとはページに書いてある内容に従ってクリックしていけば結果が出ます…と言いたいところですがいくつか注意点があります。
###5.1. まずはちゃんとQEAを選ぶ
見づらいですが3つ目のタブがQEAです。ID typeはCompound names、Data formatはSamples in columnsにしておきましょう。
図1：csvのアップロード時のオプション
###5.2. 代謝物の名称を変更する
こちらが渡した代謝物のうち、MetaboAnalystのデータベース上に完全に存在しない代謝物は解析に使われず、データベース内の名称と一致するか疑わしい代謝物はユーザー側に確認を取る形になっています。
DLが表記されてたりなかったりなど色々ありますが、大きな心をもって変更しましょう。
図2：MetaboAnalystが不安だとこんな感じで黄色の蛍光が引かれます
図3：「お前が言いたいのはこれか？」とマッチしそうなのを探してきてくれます
###5.3. サンプルの正規化とスケーリング
Quantile NormalizationとAuto-scalingを選んでおいてください。
図4：分布の正規化・スケーリングのオプション
サンプル間の分布の違いはSample Normalizationの手法に基づき補正されます。サンプル間の分布が大きく違う場合はそもそも測定がアレなのでは
代謝物は存在量のレンジが各代謝物で滅茶苦茶違うのでスケーリングしておきます。元の分布(左)だとSucroseとGlucoseが飛びぬけて存在していますが、スケーリング後(右)は分布が他の代謝物と似た感じになってるのが見て取れます。
図5：スケーリング前(左)とスケーリング後(右)の結果(一部)
###5.4. Metabolite SetはPathway-associated (KEGG)で
用意されているMetabolite Setのどれを選ぶかですが、KEGGから取得されたこれにしておきましょう。
一番上のSMPDBベースのMetabolite Setは細かすぎて逆に分かりづらいので…(あとHuman-specificっぽい？)
図6：Metabolite Set選択時のオプション
###5.5. 結果
結果はこんな感じになります。ちゃんとMethionine関連のパスウェイが一番変化したと出力していますね。
図7：MSEA結果のOverview
#6．MetaboAnalystRについて
WebのMetaboAnalystで上記の手順を行うと、右側にRのスクリプトが自動で生成されているのがわかると思います。
再現性確保や諸々のためにMetaboAnalystR5としてRパッケージが管理されています。
コードを保存しておくと、いつか解析をやり直す際に便利です。
内部的にCairo6を使っているらしく、MetaboAnalystRをローカルにインストールするには
- 必要なライブラリとCairoのインストール
- 必要なRパッケージのインストール
- MetaboAnalystRのインストール
...~~と非常にめんどくさいので、大人しくDocker7を使うことも可能です。~~Dockerfileしかなくあまりメンテされてないようなので、耐え難きを耐えながらローカル環境に整備しちゃった方がいいのかもしれません。
以下、僕が解析を行った際のコードを付記しておきます。
mSet<-InitDataObjects("conc", "msetqea", FALSE)
mSet<-Read.TextData(mSet, "gctableme.csv", "colu", "disc");
mSet<-CrossReferencing(mSet, "name");
mSet<-CreateMappingResultTable(mSet)
mSet<-PerformDetailMatch(mSet, "2-Hydroxyglutaric acid");
mSet<-GetCandidateList(mSet);
mSet<-SetCandidate(mSet, "2-Hydroxyglutaric acid", "L-2-Hydroxyglutaric acid");
mSet<-PerformDetailMatch(mSet, "Dehydroascorbic acid dimer");
mSet<-GetCandidateList(mSet);
mSet<-SetCandidate(mSet, "Dehydroascorbic acid dimer", "Dehydroascorbic acid");
mSet<-PerformDetailMatch(mSet, "D-Mannose-6-phosphate");
mSet<-GetCandidateList(mSet);
mSet<-SetCandidate(mSet, "D-Mannose-6-phosphate", "Mannose 6-phosphate");
mSet<-PerformDetailMatch(mSet, "D-Glucose-6-phosphate");
mSet<-GetCandidateList(mSet);
mSet<-SetCandidate(mSet, "D-Glucose-6-phosphate", "Beta-D-Glucose 6-phosphate");
mSet<-SanityCheckData(mSet)
mSet<-ReplaceMin(mSet);
mSet<-PreparePrenormData(mSet)
mSet<-Normalization(mSet, "QuantileNorm", "NULL", "AutoNorm", ratio=FALSE, ratioNum=20)
mSet<-PlotNormSummary(mSet, "norm_0_", "png", 72, width=NA)
mSet<-PlotSampleNormSummary(mSet, "snorm_0_", "png", 72, width=NA)
mSet<-SetMetabolomeFilter(mSet, F);
mSet<-SetCurrentMsetLib(mSet, "kegg_pathway_new", 2);
mSet<-CalculateGlobalTestScore(mSet)
mSet<-PlotQEA.Overview(mSet, "qea_0_", "net", "png", 72, width=NA)
mSet<-SaveTransformedData(mSet)
実行すると以下のファイルが得られます。
- 図5で出した正規化・スケーリング前後の分布
- norm_0_dpi72.png
- snorm_0_dpi72.png
- 入力に使用した代謝物名と、MetaboAnalyst側で認識した名称、他の公共データベース上のIDの対応表
- name_map.csv
- MSEAの実行結果(画像とcsv)
- qea_0_dpi72.png
- msea_qea_result.csv
- オリジナルのデータ、正規化・スケーリング後のデータ
- data_original.csv
- data_processed.csv
- data_normalized.csv
#7．おわりに
こんな感じで代謝物のパスウェイに基づくエンリッチメント解析手法をバーッと紹介しました。
これを読んだ人の中で現在代謝物解析をしている人、将来代謝物解析に関わる人がどれだけいるか分かりませんが、何かのきっかけになれば幸いです。
細かいところをすっ飛ばしていたり僕の説明が至らない点があったりしますので、もし何かあればお気軽にコメント下さい。