すごく更新が遅くなりました.およそ,1ヶ月ぶり…すみません…
さて,保持時間予測のお話ですが,使う武器は,
1.保持時間がわかっている化合物構造ファイル(SDF)と,保持時間を予測したい化合物のSDFファイル
2.化合物構造式と化合物プロパティの一覧表
3.PLS-Rの知識
https://sites.google.com/site/esitomonokai/jie-xi-bu-wu/tahenryo/plsniyoru-hui-gui
です.
この3つがちゃんとできていれば,SIMCA-P等,多変量解析ソフトウェアを持っていれば,もはやボタンを押すだけで保持時間予測はできると思います.一応,フリーの多変量解析ツールで,今回は紹介しようと思いますが...
管理人が脂質分析の研究で保持時間の予測をやったときは,254種類の保持時間がわかっている脂質のSDFからPaDEL-Descriptorで2000弱の変数をまず出力しました(CDKやその他の関数で出力できる2D化合物プロパティに加えて,PubChem fingerprintの値を含む).
それからまず,254の脂質すべてにおいて,「変数の標準偏差が0」の変数はすべてエクセルから削除しました.この時点で,変数の数は464種類になりました.
つまり,254脂質分子種の保持時間セットを,464種類の変数を用いてPLS回帰しましょうという問題となりました.
PLS-Rに関しては,下記リンクからダウンロードできる多変量解析用エクセルマクロを用いてやりました.
http://prime.psc.riken.jp/Metabolomics_Software/StatisticalAnalysisOnMicrosoftExcel/index.html
このエクセルマクロで,所定のフォーマット(行に変数名,列にサンプル名が来るように.クラス名は適当に全部1を入れました.)でエクセルの内容を貼り付けて,変数名の一番下の行にY変数として,保持時間情報を入力しました.
あとは,エクセルマクロのボタンに従ってX変数(254種類)とY変数(保持時間)を選択して,PLS-RのAnalysisボタンを押しただけです.
ちなみに,TransformはNoneを選択し,ScaleはAuto scaleを選択しました.
Analysisを押すと,このプログラムは「いくつまで潜在変数を計算しますか?」と聞いてくるので,管理人は20と入力し,Doneを押しました.
SIMCA-Pなら,最適な潜在変数の数を自動的に割り出してくれるのですが,このプログラムは自分で考えなければなりません.PLS-Rの最適な潜在変数の数は,Y変数が1つの場合,だいたい5前後となりますので(経験則),適当に20を入れるようにしています.
このプログラムは7-Foldのクロスバリデーションを自動的に行なってくれて(この辺はSIMCA-Pと一緒),それぞれの潜在変数の数におけるPRESS値とQ2値も出してくれます.
PRESS値とQ2値って何?という質問に対して,このコラムで答えるのはしんどいので,データの解釈だけ.
1.PRESS値は理想的には,「最適な潜在変数の数」の時に極小値を取る
2.Q2値は理想的には,「最適な潜在変数の数」を1つ超える負の値となる
上記理由に従い,潜在変数20個まで計算した中で,PRESSの極小値,もしくはQ2値が正の値から負の値となったときの潜在変数の数を選択すれば良い!ということになるのですが,実際,PRESS値が「極小値」を取ることは稀で「なだらかに減少する」ことが多く,Q2値も,「正の値だけどほとんど0じゃね?」といったことが多いです.
管理人の実際のデータですが,このとき,
というような感じでした.
PRESS値はなだらかに減少し続け,Q2値は潜在変数の数が6を超えたあたりから0.05前後と非常に小さくなっていました.
管理人は,このとき,最適な潜在変数の数を「6」として,このモデルを最適化しました.
なぜ6かといいますと,
1.実際Q2値が負となる10の1つ手前の9だと,なんかオーバーフィットしてそうな気がする.
2.Q2値が6以降,ほとんど変わっていない.
という理由から,6と決定しました.
主観が入っています.正直.
この辺をしっかり,客観的にやる研究もされていますが,僕はどちらかというと多変量解析は「ユーザー」の立場なので,
「これで作ったモデルが実際の未知データにちゃんと予測精度を保てているならそれで良い」
と思ってしまうほうなので…
次に,モデルの精度を表すRMSEEの話ですが,管理人はRMSEEは直感的ではなさすぎるので嫌いで,普通に「予測値」と「実測値」の差分から「標準偏差」を計算して,それで評価するようにしています.
このときのモデルは15分の液体クロマトグラフィーの分析系で,標準偏差が0.14分,つまり8秒前後の誤差で保持時間予測可能なモデルをトレーニングセットで作ったことになります.
次に,標準品情報が無い,保持時間情報が未知の全117,343種類の脂質に関して,同様に変数の数が464種でサンプルの数が117343のエクセルファイルを作成し,それを多変量解析マクロの「PLS-Rのトレーニングセット結果が残った状態」で,Predictionシートにその464×117343のデータを貼り付け,Predictionボタンを押して,全保持時間を出力しました.
このとき,新たに同定できた(標準品情報として入っていたものも,さらに同定の重複も含みますが)全1808種の脂質の実測保持時間と,予測保持時間の標準偏差は0.14分でしたので,オーバーフィットなくモデルが作れていたかなって思います.
以上で今回の保持時間予測のお話は終わります.
次回からはしばらく,特にテーマを決めることなく,共有したい「ちょっとしたコラム」をたてつづけに書いていこうと思います.