質量分析データから分子構造を推定するKaggleコンペに挑戦しました。
ただ、参加した時点では「そもそも何を予測するコンペなのか」がよく分かっていなかったため、最初にコンペの説明文から確認しました。
どうやら、質量分析によって得られたデータから分子構造を推定するモデルを作るコンペで、学習データに含まれている分子だけでなく、未知の分子構造についても予測することが目的のようでした。
分子構造の決定というと、1HNMRでは原子同士のつながりや周囲の化学的環境など、かなり細かな情報を得ることができます。複数の測定を組み合わせれば、立体構造について調べることもできます。
一方、今回扱うのは質量分析です。
質量分析では、イオンの質量と電荷の比である「m/z」を測定します。また、分子を壊したときに生じる断片についても、それぞれのm/zや観測強度を調べることができます。
そのため、
「元の分子がどのくらいの重さなのか」
「壊したときに、どのくらいの重さの断片ができるのか」
といった情報は得られますが、それだけから分子構造を一意に決定できるとは限りません。
今回の評価では立体異性まで区別する必要はなく、原子同士のつながりを表す二次元構造を予測できればよいため、まずは実際のデータを確認してみることにしました。
pandasのDataFrameで学習データを見ると、一部の列には通常の数値ではなく、NumPyのndarrayが格納されていました。
さらに確認すると、この配列の長さは行によって異なっていました。
その中で、ある2列のndarrayは各行で配列長が一致していました。
列名を確認し、実際に配列内部の値も見てみると、
・片方が、分子を壊したときに観測された断片のm/z
・もう片方が、それぞれのm/zに対応するスペクトル強度
という対応関係になっていることが分かりました。
例えば、
m/z = 91 → 強度 0.82
m/z = 105 → 強度 0.31
というように、ある質量付近の断片がどの程度強く観測されたかを表しています。
スペクトルの中身を確認したとき、最初はメチル基なら約15、ヒドロキシ基なら約17、エチル基なら約29といった、小さな部品の質量が多く現れるのではないかと考えていました。
ところが実際のデータを見ると、そのような小さな値はほとんど見当たりませんでした。低いものでもm/zが80〜90程度あり、メチル基などと比べてもかなり大きな値です。
ここで、自分が最初にイメージしていた「元の分子からちぎれて飛んでいった小さな部品の質量」を見ているわけではないことに気づきました。
観測されているのは、分子が切断された後も電荷を持って残った、比較的大きな断片側です。
つまり、
「何がちぎれて飛んでいったのか」
ではなく、
「ちぎられた結果、どのくらいの大きさの断片が残ったのか」
を直接見ていることになります。
そこで逆に、元の分子に近いイオンのm/zから観測された断片のm/zを引けば、
「切断によってどの程度の質量が失われたのか」
という別の情報を作れるのではないかと考えました。
これが後で、元のイオンと断片の質量差も特徴量として追加するきっかけになりました。
一方、このままでは1行ごとに配列長が異なるため、一般的な表形式の機械学習モデルにはそのまま入力できません。
そこで、m/zを列名、スペクトル強度をその列の値として、可変長の配列を固定長の表形式データへ展開することにしました。
例えば、m/z 91付近のピークであれば、
mz_91 = 0.82
という特徴量にします。
そのスペクトルにm/z 91付近のピークが存在しなければ0とします。
元のm/zは小数を含んでいますが、最初のベースラインでは細かい質量差をいったん無視し、整数単位に丸めることにしました。
炭素原子の質量がおよそ12、水素原子がおよそ1であることを考えると、まずは整数単位でも分子断片の大まかな違いを捉えられるのではないかと考えたためです。
また、1スペクトルに含まれるピーク数を確認すると、少ないものから非常に多いものまであり、多いものでは数万個のピークを含んでいました。
全てのm/zをそのまま列にすると特徴量数が膨大になるため、まずは学習データ全体で頻繁に現れるm/zだけに絞ることにしました。
こうして、
「頻出するm/zを特徴量にする」
「そのm/zが観測されたときのスペクトル強度を値にする」
という単純な表形式データを作りました。
一方、学習データ全体はかなり大きく、全ての分子構造をそのまま正解ラベルとして教師あり学習させるのも現実的ではありませんでした。
そこで、まずは出現頻度の高い100種類の分子だけを個別の正解ラベルとして残し、それ以外の分子を一旦まとめて「UNKNOWN」という1つのラベルとして扱いました。
つまり、
「このスペクトルは頻出100種類のどれか」
「それとも、現在のモデルでは分からない分子か」
という101クラスの問題に単純化しました。
最初から未知分子の構造そのものを生成するのではなく、まずはスペクトルから既知分子を識別できるか、さらに未知分子を未知として区別できるかを確認するためです。
正解として格納されている分子構造の文字列も、最初に見たときはかなり意味不明に見えました。
しかし、よく見るとSMILESという形式で、有機化学で見慣れた構造を文字列で表現していることが分かりました。
炭素鎖、芳香族環、カルボニル基、ヒドロキシ基などが、一定の文字や記号の組み合わせとして記述されています。
そこで、
「分子全体を1つの巨大な正解ラベルとして扱わなくてもよいのではないか」
と考えるようになりました。
主鎖、芳香族環、官能基などを部品として考え、それらをトークンとして扱えば、最終的に組み合わせて分子構造を再構成できる可能性があります。
分子全体の種類は非常に多い一方で、頻繁に登場する官能基や部分構造の種類はそれほど無制限ではありません。
場合によっては、現在設定している分子単位の正解ラベル数よりも少ない種類の部品に分解できる可能性があります。
このようなトークン列を扱うモデルとして代表的なのが、現在LLMでも広く使われているTransformerです。
ただし、ニューラルネットワークで構造生成まで行うと、一つの仮説を試すための学習時間も長くなります。
まだ「どの情報が有効なのか」を調べている段階で毎回大きなニューラルネットワークを学習させるのは効率が悪いため、まずはLightGBMのような比較的単純な機械学習モデルを使って実験を繰り返すことにしました。
最終的にSMILESを生成する段階ではTransformerなどを検討するとしても、その前に、
「質量そのものがどの程度効くのか」
「分子が壊れたときの情報はどの程度効くのか」
「測定条件に依存した情報を拾っていないか」
を切り分けておきたいと考えました。
最初に、元の分子に近いイオンの質量と、頻出する断片のm/zを100種類使って学習したところ、Accuracyは約87.1%でした。
さらに、測定時に分子へ何が付加しているかを表すadduct情報を追加すると、Accuracyは約90.2%まで上昇しました。
同じ分子でも、水素イオンが付加した[M+H]+とナトリウムイオンが付加した[M+Na]+では観測される質量が変わるため、これは化学的にも納得できる結果でした。
そこで、adductを利用して、観測された質量から元の分子そのものの質量に近い値を計算し、推定中性分子量として追加しました。
さらに、スペクトルを確認した際に、
「観測されているのは、ちぎれて飛んでいった小さな部品ではなく、ちぎられたあとに残った側の断片なのではないか」
と気づいたことから、
「元の分子に近いイオンの質量 − 壊れてできた断片の質量」
も計算してみることにしました。
これなら、直接観測されている断片側の質量から、逆に「元の分子からどのくらいの質量が失われたのか」を表現できます。
例えば、元のイオンと断片の質量差が約18であれば、水分子程度の質量が失われた可能性などを考えることができます。
そこで、
・元の分子に近いイオンの質量
・測定時に付加したイオンの種類
・推定した中性分子量
・頻出する断片のm/z 500種類
・元のイオンと断片の質量差として頻繁に現れる値500種類
を合わせた約1000個の特徴量を作りました。
この特徴量を使ったモデルではAccuracyが約95.6%まで上昇しました。
ただし、学習を続けると途中から誤差が急激に悪化する場合があり、モデルが複雑になりすぎている可能性がありました。
そこでLightGBMの木を小さくし、L1・L2正則化やデータ・特徴量のサンプリングを加えて過学習を抑えました。
その結果、現在の基準モデルでは、
Accuracy：約97.2%
MRR@25：約0.985
UNKNOWN判定AUC：約0.997
Logloss：約0.094
となっています。
また、m/zを整数単位ではなく0.1単位まで細かく分ける実験も行いました。
炭素原子の質量がおよそ12、水素原子がおよそ1であることを考えると、0.1単位まで細かくしたからといって大きく性能が上がるとは考えていませんでした。
ただし、実際の原子や同位体の質量は完全な整数ではなく、細かな質量差にも何らかの情報が含まれている可能性があります。
そのため、「整数に丸めることで使える情報まで失っていないか」を確認する目的で、0.1単位に細かくした特徴量も一度試してみました。
結果として、整数単位のm/zを使った場合のAccuracyが約90.2%だったのに対し、0.1単位では約86.9%となり、今回のモデルでは性能が低下しました。
少なくとも現在のLightGBMによる分類では、細かな質量差をそのまま多数の特徴量へ分けるより、整数単位でまとめた方が扱いやすい結果となりました。
現在は、約1000個の特徴量をそのまま使うだけでなく、
・推定中性分子量を削除する
・付加イオン情報を削除する
・断片m/zを削除する
・元のイオンと断片の質量差を削除する
・使用するm/zや質量差を100、200、500種類に変える
といった実験を行い、どの情報が本当に性能へ寄与しているのかを一つずつ確認しています。
また、通常のランダムな学習・検証分割では、同じ分子の別スペクトルが学習側と検証側の両方に入る可能性があります。
そこで、UNKNOWNに分類した分子については、分子構造を表すInChIKey14単位で丸ごと学習側と検証側に分離し、学習時に一度も見ていない分子に対してもUNKNOWNと判断できるかを確認する実験も行っています。
この実験では、正則化前のモデルでもUNKNOWN判定AUCが約0.979となりました。
もちろん、これは「未知分子のSMILESそのものを正しく生成できた」という意味ではありません。
現時点で確認しているのは、
「学習した100種類のどれでもなさそうだ」
と判定できるところまでです。
今後は、現在行っている特徴量の比較を終えた後、分子全体を1つのラベルとして分類する方法から、官能基や部分構造などを予測する問題へ少しずつ移し、最終的にはそれらを組み合わせて未知分子の構造を生成する方法を試したいと考えています。
今回取り組んでいて面白いのは、単にモデルを複雑にすることではなく、
「測定された数値が化学的に何を表しているのか」
「どのように変換するとモデルが情報を捉えやすくなるのか」
「実際に削除したとき性能が下がるのか」
を一つずつ実験で確認できるところです。
現在も特徴量の削除実験や、使用するm/z・質量差の数を変える比較、未知分子をより厳密に分離した検証などを続けています。
今後、実験結果がまとまり次第、このストーリーにも順次追加していく予定です。
