特許 7440870 ２つの変数により定まるデータの解析方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-02-20

(45)【発行日】2024-02-29

(54)【発明の名称】２つの変数により定まるデータの解析方法

(51)【国際特許分類】

   G01N 27/62 20210101AFI20240221BHJP

   G01N 21/35 20140101ALI20240221BHJP

   G01N 30/72 20060101ALI20240221BHJP

【ＦＩ】

G01N27/62 D

G01N21/35

G01N30/72 C

【請求項の数】 9

(21)【出願番号】P 2021009134

(22)【出願日】2021-01-22

(65)【公開番号】P2022113039

(43)【公開日】2022-08-03

【審査請求日】2023-05-09

(73)【特許権者】

【識別番号】000001993

【氏名又は名称】株式会社島津製作所

(73)【特許権者】

【識別番号】504132272

【氏名又は名称】国立大学法人京都大学

(74)【代理人】

【識別番号】110001069

【氏名又は名称】弁理士法人京都国際特許事務所

(72)【発明者】

【氏名】小澤 弘明

(72)【発明者】

【氏名】藤田 雄一郎

(72)【発明者】

【氏名】石濱 泰

(72)【発明者】

【氏名】吉井 和佳

【審査官】吉田 将志

(56)【参考文献】

【文献】特開２０１８－００４２５２（ＪＰ，Ａ）

【文献】特開２０１６－０５７０６５（ＪＰ，Ａ）

【文献】特開２０１６－１３０７２７（ＪＰ，Ａ）

【文献】国際公開第２０１７／０７７６１８（ＷＯ，Ａ１）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｎ ２７／６０ － Ｇ０１Ｎ ２７／７０

Ｇ０１Ｎ ２１／００ － Ｇ０１Ｎ ２１／９５８

Ｇ０１Ｎ ３０／００ － Ｇ０１Ｎ ３０／９６

(57)【特許請求の範囲】

【請求項1】

分析装置から得られた、第１の変数の値が互いに異なるN行のデータと第２の変数の値が互いに異なるM列のデータから成るN行M列のデータ行列Ｘに基づいて、N行K列の第１行列ＳとK行M列の第２行列Ｐの積ＳＰで該データ行列Ｘを近似するように、該第１行列Ｓ及び該第２行列Ｐを求める方法であって、
正則化パラメータの候補λr（rは1からr_maxまでの自然数）を複数個、及び解のスパース性を誘導する正則化関数R(Ｓ, Ｐ)を1個用意する正則化パラメータ・正則化関数準備工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記データ行列Ｘと前記積ＳＰの相違度を示す距離関数D(Ｘ｜ＳＰ) と、正則化パラメータ候補λrと正則化関数R(Ｓ, Ｐ)の積λrR(Ｓ, Ｐ)との和である損失関数 L(Ｓ, Ｐ) =D(Ｘ｜ＳＰ)+λrR(Ｓ, Ｐ)の値が最小となるように求められる行列Ｓrt及び行列Ｐrtを第１行列Ｓの候補Ｓr及び第２行列Ｐの候補Ｐrとして決定するように最適化問題を解く行列候補決定工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記データ行列Ｘの行列要素X_nm及びそれに対応する前記第１行列候補Ｓrと前記第２行列候補Ｐrの積ＳrＰrの行列要素(SrPr)_nmの組み合わせ毎に、前記距離関数D(X_nm｜(SP)_nm)に対応する確率分布P_nmを共通の確率分布 Pcommonへ変数変換する関数F_nmを用いて、共通の確率分布Pcommonへと変数変換された値である変換値y_nm=F_nm(X_nm｜(SrPr)_nm)を求める確率分布変換工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記変換値y_nmと確率分布Pcommonの累積分布関数との適合度を求める適合度算出工程と、
前記複数個の正則化パラメータの候補λrのうち、前記適合度が最大値となるもの、又は前記適合度が所定の閾値以上であって前記第１行列候補Ｓr及び前記第２行列候補Ｐrにおける前記λrの値が最大であるものについて決定した前記第１行列候補Ｓr及び前記第２行列候補Ｐrを、第１行列Ｓ及び第２行列Ｐとして決定する行列決定工程と
を有するデータ解析方法。

【請求項2】

前記行列候補決定工程において、前記行列Ｓrt及び前記行列Ｐrtを第１行列候補Ｓr及び第２行列候補Ｐrと決定する代わりに、該行列Ｓrt及び該行列Ｐrtを初期値として、正則化項を含まない第２の損失関数 L₂(Ｓ, Ｐ)= D(Ｘ|ＳＰ) を最小化するように行列Ｓrt2及び行列Ｐrt2を求め、該行列Ｓrt2及び該行列Ｐrt2を第１行列候補Ｓr及び第２行列候補Ｐrと決定する、請求項１に記載のデータ解析方法。

【請求項3】

前記正則化関数Ｒ(Ｓ, Ｐ)が、L1ノルム、若しくはL1ノルムとL2ノルムの線形結合、又は、Ｐの各列の総和が1以下になるように解の制約を加えたうえでＳに対してトレースノルム、det|Ｓ^TＳ|、logdet|Ｓ^TＳ+δI|（但し、Iは単位行列、δは正則化関数を制御するハイパーパラメータ）のいずれかを適用したものである、請求項１又は２に記載のデータ解析方法。

【請求項4】

前記累積分布関数が、ガウス分布を表す関数から求められる累積分布関数、ポアソン分布を表す関数から求められる累積分布関数、指数分布を表す関数から求められる累積分布関数、Tweedie分布を表す関数から求められる累積分布関数のいずれかである、請求項１～３のいずれか１項に記載のデータ解析方法。

【請求項5】

前記適合度がコルモゴロフ－スミルノフ検定、クラーメル－フォン ミーゼス検定、及びアンダーソン－ダーリング検定のうちのいずれかを用いて求められる、請求項１～４のいずれか１項に記載のデータ解析方法。

【請求項6】

前記行列要素X_nm毎のノイズの分散σ_nm ²が既知である場合において、前記確率分布Pcommonを標準正規分布とし、前記累積分布関数F_nm(X_nm｜(SrPr)_nm)を (X_nm - (SrPr)_nm)/σ_nmとし、前記変換値y_nm が平均0であると想定した場合の不偏分散

【数8】

を用いて表される-|σ_y ² - 1|を前記適合度として用いる、請求項１～４のいずれか１項に記載のデータ解析方法。

【請求項7】

前記分析装置がクロマトグラフ質量分析装置であり、前記第１の変数及び前記第２の変数のいずれか一方がm/zであって他方が保持時間である、請求項１～６のいずれか１項に記載のデータ解析方法。

【請求項8】

前記分析装置がイメージング質量分析装置であり、前記第１の変数及び前記第２の変数のいずれか一方がm/zであって他方が位置を表す数値である、請求項１～６のいずれか１項に記載のデータ解析方法。

【請求項9】

前記分析装置がイメージングフーリエ変換赤外分光光度計であり、前記第１の変数及び前記第２の変数のいずれか一方が吸収波長であって他方が位置を表す数値である、請求項１～６のいずれか１項に記載のデータ解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、測定対象の成分に依存する第１の変数とその他の第２の変数という2つの変数により定まるデータを解析する方法に関する。例えば、クロマトグラフ質量分析装置で得られる保持時間毎の各m/z（質量電荷比）の検出強度のデータでは、m/z（質量電荷比）が第１の変数に該当し、保持時間が第２の変数に該当する。イメージング質量分析装置で得られる位置毎のm/z強度のデータでは、m/zが第１の変数に該当し、位置が第２の変数に該当する。イメージングFT-IR（フーリエ変換赤外分光光度計）で得られる位置毎の吸収スペクトルのデータでは、吸収波長が第１の変数に該当し、位置が第２の変数に該当する。

【背景技術】

【0002】

上記のように2つの変数により定まるデータは、第１の変数の値が互いに異なるN行のデータと、第２の変数の値が互いに異なるM列のデータから成るN行M列のデータ行列で表される。ここでN, Mはそれぞれ自然数である。例えば、液体クロマトグラフ質量分析装置（ＬＣ／ＭＳ）やガスクロマトグラフ質量分析装置（ＧＣ／ＭＳ）といったクロマトグラフ質量分析装置では、或る保持時間においてN種類のm/zにおける検出強度を示すN個のデータから成るマススペクトルが取得され、それらのデータがデータ行列中の或る1個の列におけるN個のデータとなる。そして、異なるM個の保持時間においてそれぞれマススペクトルが取得されることにより、M個の列においてそれぞれN個のデータが得られる。その結果、N行M列のデータ行列が得られる。

【0003】

クロマトグラフ質量分析装置により測定される試料には1又は複数種類の成分が含まれている。なお、ここで試料に含まれている「成分」とは、クロマトグラムやマススペクトルに所定以上の高さ（又は大きさ）のピークを形成する有意な成分のことを言う。試料にK種類（Kは自然数）の成分が含まれているとすると、保持時間毎のマススペクトルやm/z毎のマスクロマトグラムは、それらK種類の成分からそれぞれ得られるマススペクトルやマスクロマトグラムを重ね合わせたものになる。そのように重ね合わされたマススペクトルを成分毎のK種類のマススペクトルに分離したり、重ね合わされたマスクロマトグラムを成分毎のK種類のクロマトグラムに分離することができれば、成分毎の情報を得ることができる。そのために、行列分解と呼ばれる手法を用いることができる（例えば特許文献１）。

【0004】

行列分解では、前述したN行M列のデータ行列Ｘを、N行K列のスペクトル行列ＳとK行M列のプロファイル行列Ｐの積ＳＰで表すことができるように、スペクトル行列Ｓ及びプロファイル行列Ｐを求める。ここでKは上述のように試料に含まれる成分の種類の数に相当するが、行列分解の手法においては「因子数」と呼ばれる。スペクトル行列Ｓの各列はK種類の成分毎のマススペクトルに対応し、プロファイル行列Ｐの各行はK種類の成分毎のクロマトグラムに対応する。

【0005】

これらスペクトル行列Ｓ及びプロファイル行列Ｐは、通常、解析的に求めることはできない。そこで、コンピュータを用いて、複数のスペクトル行列Ｓ及びプロファイル行列Ｐの候補に対して、各スペクトル行列Ｓの候補及び各プロファイル行列Ｐの候補の積ＳＰについてそれぞれ、行列要素毎に目的とするデータ行列Ｘで対応する行列要素との誤差を所定の関数 D(x|y)（「距離関数」または「ダイバージェンス」と呼ばれる）により定義したうえで、それら行列要素毎の誤差の総和を求め、その値が最小になるようにスペクトル行列Ｓ及びプロファイル行列Ｐを求める、という手法が用いられている。距離関数D(x|y)には(x-y)²で表されるユークリッド距離、xlog(x/y)-(x-y)で表される一般化ＫＬ（Kullback-Leibler）ダイバージェンス、(x/y)-log(x/y)-1で表されるＩＳ（Itakura-Saito）ダイバージェンス、それらを一般化したβダイバージェンス

【数1】

等が用いられ、この関数の行列要素毎の総和 D(Ｘ|ＳＰ) = Σ_n,m D(X_nm|(SP)_nm)がデータ行列Ｘと積ＳＰの相違度を示すこととなる。ここで(SP)_nmは積ＳＰから成る行列の行列要素である。

【0006】

一方で、確率分布のパラメータの推定問題では、最尤推定によってパラメータの推定値を得ることがよく行われている。すなわち、パラメータθ(θは複数の次元を持つベクトルであってもよい)をもつ確率分布において、データxに対するパラメータθの尤(もっと)もらしさを表す尤度関数 p(x|θ)が最大となるようにθを求める。ここで、尤度の最大化により得られるパラメータθの推定値と、負の対数尤度関数 -log p(x|θ) を最小化して得られるパラメータθの推定値は一致するため、より計算の行いやすい負の対数尤度の最小化がよく用いられる。

【0007】

さて、行列分解を最尤推定問題として解釈する場合、上記の距離関数 D(Ｘ|ＳＰ) を最小化するように行列Ｓ, Ｐを求めることは、行列の各要素毎に距離関数と対応する確率分布と尤度関数p(X_nm|(SP)_nm) を想定することと一致する。すなわち、確率分布として、(i)ガウス分布、(ii)ポアソン分布、(iii)指数分布、(iv)Tweedie分布、を挙げると、それぞれの負の対数尤度関数が、 (i)ユークリッド距離、(ii)一般化ＫＬダイバージェンス、(iii)ＩＳダイバージェンス、(iv)βダイバージェンス、に対応する（例えば非特許文献１）。なお、ここでは4つの例について確率分布と距離関数の対応関係を示したが、一般に、任意の確率分布について、その負の対数尤度関数を距離関数とすれば、同様の推定が可能となる。

【0008】

また、スパースモデリングの分野では、パラメータの推定値がスパース（sparse：日本語で「まばらな」との意味）となる（行列では0である行列要素が多くなる）ように、解のスパース性を誘導する正則化項としてL1ノルム等を距離関数に加えた損失関数を定義し、その損失関数を最小化することで、不要なパラメータが0と推定されるようにするLASSO（Least Absolute Shrinkage and Selection Operator）等の手法が用いられる。

【先行技術文献】

【特許文献】

【0009】

【文献】米国特許公開公報第2009/0121125号

【非特許文献】

【0010】

【文献】吉井和佳 他3名、「NMF vs PLCA: 多重音生成過程のための無限因子モデルと無限混合モデル」、情報処理学会研究報告 MUS（音楽情報科学）、vol. 2016-MUS-112、No. 21、pp. 1-10、2016年8月1日

【発明の概要】

【発明が解決しようとする課題】

【0011】

データ行列Ｘに対して行列分解を行う際には、因子数Kを定めなければならない。その際、試料に含まれる成分の種類の数が既知である場合にはその数を因子数Kとすればよいが、多くの場合にはそのような数は未知である。このように試料に含まれる成分の種類の数が未知であると、因子数Kを適切に定めることが困難であるため、データ行列Ｘに基づいてスペクトル行列Ｓ及びプロファイル行列Ｐを適切に求めることも難しい。また、スパースモデリングのように、解のスパース性を誘導する正則化項を損失関数に追加する場合でも、正則化パラメータの値を適切に決定することが難しい。

【0012】

本発明が解決しようとする課題は、試料に含まれる成分の種類の数が未知であっても、該試料に対する測定により得られたデータ行列Ｘに基づいて、適切な、すなわち、該成分の数に近い因子数Kを有するスペクトル行列Ｓ及びプロファイル行列Ｐを求めることができるデータ解析方法を提供することである。

【課題を解決するための手段】

【0013】

上記課題を解決するために成された本発明に係る2つの変数に依存するデータの解析方法は、分析装置から得られた、第１の変数の値が互いに異なるN行のデータと第２の変数の値が互いに異なるM列のデータから成るN行M列のデータ行列Ｘに基づいて、N行K列の第１行列ＳとK行M列の第２行列Ｐの積ＳＰで該データ行列Ｘを近似するように、該第１行列Ｓ及び該第２行列Ｐを求める方法であって、
正則化パラメータの候補λr（rは1からr_maxまでの自然数）を複数個、及び解のスパース性を誘導する正則化関数R(Ｓ, Ｐ)を1個用意する正則化パラメータ・正則化関数準備工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記データ行列Ｘと前記積ＳＰの相違度を示す距離関数D(Ｘ｜ＳＰ) と、正則化パラメータ候補λrと正則化関数R(Ｓ, Ｐ)の積λrR(Ｓ, Ｐ)との和である損失関数 L(Ｓ, Ｐ) =D(Ｘ｜ＳＰ)+λrR(Ｓ, Ｐ)の値が最小となるように求められる行列Ｓrt及び行列Ｐrtを第１行列Ｓの候補Ｓr及び第２行列Ｐの候補Ｐrとして決定するように最適化問題を解く行列候補決定工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記データ行列Ｘの行列要素Xnm及びそれに対応する前記第１行列候補Ｓrと前記第２行列候補Ｐrの積ＳrＰrの行列要素(SrPr)_nmの組み合わせ毎に、前記距離関数D(X_nm｜(SP)_nm)に対応する確率分布P_nmを共通の確率分布 Pcommonへ変数変換する関数F_nmを用いて、共通の確率分布Pcommonへと変数変換された変換値y_nm=F_nm(X_nm｜(SrPr)_nm)を求める確率分布変換工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記変換値y _nm と確率分布Pcommonの累積分布関数との適合度を求める適合度算出工程と、
前記複数個の正則化パラメータの候補λrのうち、前記適合度が最大値となるもの、又は前記適合度が所定の閾値以上であって前記第１行列候補Ｓr及び前記第２行列候補Ｐrにおける前記λrの値が最大であるものについて決定した前記第１行列候補Ｓr及び前記第２行列候補Ｐrを、第１行列Ｓ及び第２行列Ｐとして決定する行列決定工程と

【0014】

なお、ここでは便宜上、第１行列には前記スペクトル行列と同じ記号「Ｓ」を付し、第２行列には前記プロファイル行列と同じ記号「Ｐ」を付したが、記号の通りに第１行列をスペクトル行列とし、第２行列をプロファイル行列としてもよいし、それとは反対に、第１行列をプロファイル行列とし、第２行列をスペクトル行列としてもよい。後者の場合、前述の定義とは異なり、スペクトル行列がK行M列の行列として、プロファイル行列がN行K列の行列として定義される。

【0015】

この方法によれば、スパース性を誘導する正則化の効果により、試料に含まれる成分の種類の数が未知であっても不要な因子である列が全て0になった第１行列候補Ｓｒ、及び、行の値がすべて0になった第２行列候補Ｐｒを用意することができる。一方、正則化パラメータの値が小さすぎると、0以外の値をとるSの列数(すなわち推定した因子数K)が多くなりすぎ、オーバーフィット(過剰適合。複雑なモデルを用いることにより、不合理で誤ったモデルであってもデータとの一致度が高くなってしまうこと。「過学習」とも言う。）が生じる。一方で、正則化パラメータの値が大きすぎると、推定した因子数Kが少なくなりすぎ、データにうまく適合できないアンダーフィット(過小適合、過小学習)が生じる。

【0016】

本発明では、行列候補決定工程において複数の正則化パラメータの候補λrについてそれぞれ1組ずつ決定した候補（Ｓr, Ｐr）から、確率分布変換工程から行列決定工程までの工程において理想分布との適合度（詳細は後述）を用いて絞り込みを行うことにより1組の第１行列Ｓ及び第２行列Ｐ（並びにそれらに対応する1つの因子数K）が得られる。これにより、過剰適合・過小適合が生じることが防止される。すなわち、因子数Kを実際の成分の数よりも多く設定してしまうことが防止され、適切な因子数Kを有する第１行列Ｓ及び第２行列Ｐを求めることができる。

【0017】

解のスパース性を誘導する正則化関数として、一般に「L1ノルム」、「L1ノルムとL2ノルムの線形結合」、「ボリューム制約」と呼ばれるものが知られている。本発明における正則化関数Ｒ(Ｓ, Ｐ)において、L1ノルムとL2ノルムの線形結合はＲ(Ｓ, Ｐ)=α(|Ｓ|₁+|Ｐ|₁)+(1-α)(|Ｓ|₂ ²+|Ｐ|₂ ²)で表され（但し、|Ｓ|₁=Σ_{i, j}Ｓ_{i, j}、|Ｐ|₁=Σ_{i, j}Ｐ_{i, j}、|Ｓ|₂=(Σ_{i, j}Ｓ_{i, j} ²)^1/2、|Ｐ|₂=(Σ_{i, j}Ｐ_{i, j} ²)^1/2、αは0から1の間の定数）で表され、L1ノルムはα=1の場合に該当する。「ボリューム制約」は、Ｐの各列の総和が1以下になるように解の制約を加えたうえで、Ｓに対しては、トレースノルム（シャッテン1ノルム）、det|Ｓ^TＳ|、logdet|Ｓ^TＳ+δI|（但し、Iは単位行列、δは正則化関数を制御するハイパーパラメータ）のいずれかを適用したものである。

【0018】

確率分布の累積分布関数（ここでは一般化した変数xを用いてf(x)と標記する）は、それが前提とするデータの確率分布を示す確率密度関数p(x)を用いて

【数2】

により求められる関数である。確率密度関数p(x)には、例えば、上述した(i)ガウス分布、(ii)ポアソン分布、(iii)指数分布、(iv)Tweedie分布、の各関数が挙げられる。ここで例示した4種の確率密度関数p(x)に対応する累積分布関数f(x)はそれぞれ、上述した(i)ユークリッド距離、(ii)一般化ＫＬダイバージェンス、(iii)ＩＳダイバージェンス、(iv)βダイバージェンスに、それぞれ対応している。

【0019】

確率分布変換工程で求められる変換値y_nm=F_nm(X_nm｜(SrPr)_nm)は、データの各行列要素X_nm及び積ＳrＰrの各行列要素(ＳrＰr)_nmに対応する確率分布Pnm を、共通の確率分布Pcommonへと変数変換するように求められる。例えば、一般化KLダイバージェンスを距離関数としている場合、対応する確率分布Pnmはポアソン分布となる。この時、任意の確率変数は、累積分布関数による変数変換すると標準一様分布に従うこと(確率積分変換: Probability integral transform)を利用して、F_nmを累積分布関数とすることで、Pcommon を標準一様分布とすることができる。標準一様分布に変数変換された変換値y_nmがなす経験分布関数は、行列分解の妥当性が高いほど、想定した確率分布Pcommonの累積分布関数との適合度(goodness of fit)が高くなる。そこで適合度算出工程では、複数個の正則化パラメータの候補λrのそれぞれについて、変換値y_nmと確率分布Pcommonの累積分布関数との適合度を求める。ここで適合度は、コルモゴロフ-スミルノフ統計量等の既知の適合度を利用できる。また、2つ目の具体例として、データの行列要素X_nm毎のノイズの分散σ_nm ²が既知である場合に、確率分布Pcommonを標準正規分布とし、累積分布関数F_nm(X_nm｜(SrPr)_nm)を (X_nm - (SrPr)_nm)/σ_nmとし、適合度として、変換値y_nm が平均0であると想定した場合の不偏分散

【数3】

が1に近いかどうかを表す -|σ_y ² - 1| を用いても良い。

【0020】

このように適合度を求めたうえで、行列決定工程において、(i)適合度が最大値となるもの、又は(ii)適合度が所定の閾値以上であって正則化パラメータが最大のもののいずれかについて決定した第１行列候補Ｓr及び第２行列候補Ｐrを、第１行列Ｓ及び第２行列Ｐとして決定する。ここで(ii)は、適合度が或る程度高い（所定の閾値以上である）範囲内で、できるだけ単純な（Kが小さい）モデルを用いた方が正確であると考えられることに基づいている。

【発明の効果】

【0021】

本発明により、試料に含まれる成分の種類の数が未知であっても、該試料に対する測定により得られたデータ行列Ｘに基づいて、適切な、すなわち、該成分の数に近い因子数Kを有する、スペクトル行列（第１行列又は第２行列）及びプロファイル行列（第２行列又は第１行列）を求めることができる。

【図面の簡単な説明】

【0022】

【図1】本発明の一実施形態に係るデータ解析方法が実行されるクロマトグラフ質量分析装置の一実施形態を示す概略構成図。

【図2】図１に示すクロマトグラフ質量分析装置の動作を示すフローチャート。

【図3】３次元データ、並びにマススペクトル及びクロマトグラムのデータの一例を用いて、データ行列、並びにスペクトル行列及びプロファイル行列を概念的に示す図。

【図4】図１に示すクロマトグラフ質量分析装置の動作の一部である、本実施形態のデータ解析方法の操作の詳細を示すフローチャート。

【図5】本実施形態のデータ解析方法で実行する行列分解の計算の例で用いたデータ行列を与える測定データをクロマトグラム及びマススペクトルで示す図。

【図6】データ行列を与える測定データを複数のマスクロマトグラムの重ね書きで示す図。

【図7】行列分解を実行する際に得られた、λrが小さ過ぎることにより正則化が不十分であるマススペクトル及びクロマトグラムの計算結果の例を示す図。

【図8】行列分解を実行する際に得られた、λrが大き過ぎることにより正則化が不十分であるマススペクトル及びクロマトグラムの計算結果の例を示す図。

【図9】最適なλrを用いた適切な正則化を伴う行列分解により得られた、マススペクトル及びクロマトグラムの計算結果の例を示す図。

【発明を実施するための形態】

【0023】

図１～図９を用いて、本発明に係るデータ解析方法の一実施形態、及び該データ解析方法が実行されるクロマトグラフ質量分析装置を説明する。

【0024】

(1) 本実施形態のデータ解析方法が実行されるクロマトグラフ質量分析装置の構成
図１に、本実施形態のデータ解析方法が実行される液体クロマトグラフ／イオントラップ飛行時間型質量分析装置（ＬＣ／ＩＴ－ＴＯＦＭＳ）１の要部の構成を示す。このＬＣ／ＩＴ－ＴＯＦＭＳ１は、大別して、液体クロマトグラフ（ＬＣ）部１０と、質量分析（ＭＳ）部２０と、データ処理部４０と、分析制御部５０とを有する。

【0025】

ＬＣ部１０は、移動相容器１１と、送液ポンプ１２と、インジェクタ１３と、カラム１４とを有する。移動相容器１１には移動相を貯留するものである。送液ポンプ１２は移動相容器１１内の移動相を吸引して一定流量でインジェクタ１３へ送給するものである。インジェクタ１３はオートサンプラを備え、予め用意された試料を自動的に選択し、所定のタイミングで所定料の試料を移動相中に注入するものである。インジェクタ１３により移動相中に試料が注入されると、試料は移動相の流れに乗ってカラム１４へ導入される。カラム１４を通過する間に試料中の各種成分は分離され、時間的にずれてカラム１４の出口から溶出し、ＭＳ部２０に導入される。

【0026】

ＭＳ部２０は、大気圧雰囲気に維持されるイオン化室２１と、ターボ分子ポンプ（図示せず）により真空排気され、高真空雰囲気に維持される分析室２９とを有し、その間には、段階的に真空度が高くされた第１段中間真空室２４及び第２段中間真空室２７が配設されている。イオン化室２１と第１段中間真空室２４とは細径の脱溶媒管２３を介して連通しており、第１段中間真空室２４と第２段中間真空室２７とは、円錐形状のスキマー２６の頂部に穿設された小径のオリフィスを介して連通している。第１段中間真空室２４内には第１イオンガイド２５が、第２段中間真空室２７内には第２イオンガイド２８が、それぞれ配置されている。

【0027】

イオン化室２１にはイオン源としてのＥＳＩノズル２２が配置されている。ＥＳＩノズル２２は、ＬＣ部１０から供給された試料成分を含む溶出液を、高圧電源（図示せず）から印加されている直流高電圧により帯電させた微小液滴としてイオン化室２１内に噴霧するものである。この帯電液滴は大気由来のガス分子と衝突してさらに微細な液滴に粉砕され、速やかに乾燥して（脱溶媒化されて）試料分子が気化する。この試料分子はイオン蒸発反応を生じてイオン化される。発生したイオンを含む微小液滴は、イオン化室２１と第１段中間真空室２４の差圧によって脱溶媒管２３内に引き込まれ、脱溶媒管２３内を通る間に一層脱溶媒化が進行してイオンが発生する。なお、試料分子をイオン化する方法は、ここで述べたエレクトロスプレイイオン化法（ＥＳＩ）には限られず、例えば大気圧化学イオン化法（ＡＰＣＩ）や大気圧光イオン化法（ＡＰＰＩ）を用いることができる。

【0028】

脱溶媒管２３を通過したイオンは、第１イオンガイド２５及び第２イオンガイド２８で収束されつつ第１段中間真空室２４及び第２段中間真空室２７を通過し、分析室２９へ送られる。

【0029】

分析室２９内には、イオントラップ３０、質量分離器としての飛行時間型質量分離器（ＴＯＦ）３１、及びイオン検出器３３が設けられている。

【0030】

イオントラップ３０では、電源（図示せず）より各電極に印加される高周波電圧により形成される四重極電場によって、イオンが一旦捕捉・蓄積される。イオントラップ３０の内部に蓄積された各種イオンは、所定のタイミングで一斉に運動エネルギーを付与され、イオントラップ３０からＴＯＦ３１に向けて放出される。

【0031】

また、イオントラップ３０は図示するように、例えばアルゴン等の衝突誘起解離（ＣＩＤ）ガスを供給できる構成となっており、イオントラップ３０に蓄積したイオンをＣＩＤにより開裂させてプロダクトイオンを生成させることができる。ＭＳ²分析を行う際には、まずイオントラップ３０に各種イオンを蓄積した後に、それらイオンの中で特定のm/zを有するイオンのみをプリカーサイオンとして選択的に残すように電極に印加する電圧を制御し、それからＣＩＤガスをイオントラップ３０に導入してプリカーサイオンの開裂を促進させる。そうして生成されたプロダクトイオンが所定のタイミングで一斉にイオントラップ３０からＴＯＦ３１に向けて放出される。

【0032】

ＴＯＦ３１は直流電源（図示せず）から直流電圧が印加されるリフレクトロン電極３２を備え、これにより形成される直流電場の作用によってイオンが折り返され、イオン検出器３３に到達する。一斉にイオントラップ３０から出射されたイオンはm/zの小さいイオンほど速く飛行し、m/zに応じた時間差を以てイオン検出器３３に到達する。イオン検出器３３は到達したイオン数に応じた電流を検出信号として出力する。

【0033】

イオン検出器３３には、検出信号をデジタル値に変換するＡ／Ｄ変換器３４が接続されている。Ａ／Ｄ変換器３４で変換された検出信号はデータ処理部４０に入力される。

【0034】

データ処理部４０は、データ行列作成部４１と、行列分解実行部４２と、m/z検出部４３と、保持時間検出部４４と、ＭＳ2分析実行条件決定部４５とを有する。行列分解実行部４２はその内部に正則化パラメータ・正則化関数準備部４２１と、行列候補決定部４２２と、確率分布変換部４２３と、適合度算出部４２４と、行列決定部４２５とを有する。これら各部の詳細は後述する。また、データ処理部４０は記憶部６１に接続されている。

【0035】

分析制御部５０はＬＣ／ＭＳ分析及びＬＣ／ＭＳ2分析を実行するためにＬＣ部１０及びＭＳ部２０の各部の動作を制御するものであって、ＬＣ／ＭＳ分析実行部５１及びＬＣ／ＭＳ2分析実行部５２とを有する。

【0036】

これらデータ処理部４０及び分析制御部５０は、所定の制御／処理ソフトウエアを搭載したパーソナルコンピュータ（ＰＣ）により具現化される。また、記憶部６１は当該ＰＣが有するハードディスクやソリッドステートドライブ等の記憶装置により具現化される。当該ＰＣはさらに、表示部（ディスプレイ）６２、及びキーボード、マウス、タッチパネル等のデバイスから成る操作部６３を有する。

【0037】

(2) 本実施形態のデータ解析方法を含むＬＣ／ＩＴ－ＴＯＦＭＳ１の動作
次に、図２及び図３を用いて、本実施形態のデータ解析方法を含むＬＣ／ＩＴ－ＴＯＦＭＳ１の動作を説明する。併せて、データ処理部４０が有する各構成要素の機能を説明する。

【0038】

まず、操作者が操作部６３を用いて所定の測定開始操作を行うと、ＬＣ／ＩＴ－ＴＯＦＭＳ１は、ＬＣ／ＭＳ分析実行部５１がＬＣ／ＩＴ－ＴＯＦＭＳ１の各部を制御することにより、以下のように目的試料に対するＬＣ／ＭＳ分析を実行する（ステップ１）。これにより、インジェクタ１３から移動相中に注入された目的試料がカラム１４に送り込まれ、カラム１４からの溶出液はＭＳ部２０に導入されて繰り返し質量分析される。ＭＳ部２０のイオン検出器３３から出力される検出信号は、Ａ／Ｄ変換器３４でデジタル値に変換された後に、データ処理部４０のデータ行列作成部４１に入力される。

【0039】

データ行列作成部４１では、イオントラップ３０からの1回のイオン出射によって所定のm/z範囲内でm/z毎に得られるN個の信号を、N行M列のデータ行列Ｘ中の1列に含まれるN個の行列要素の値として取得する。そして、時間差をもってなされるM回のイオン出射の各々において同様にデータを取得し、それらのデータに基づいて、下記のようなN行M列のデータ行列Ｘを作成する（ステップ２）。

【数4】

【0040】

ここで、データ行列Ｘの各行列要素X_nm（nは1～Nのいずれか、mは1～Mのいずれか）は、m回目のイオン出射の際（保持時間に対応）に検出される、前記m/z範囲内でのn番目のm/zの強度を示す。データ行列Ｘの各行列要素X_nmは0又は正の値（非負値）を有する。

【0041】

次に、行列分解実行部４２は、後述（「(3)本実施形態のデータ解析方法（行列分解の操作）の詳細」）の方法により行列分解を実行し、データ行列ＸをN行K列のスペクトル行列ＳとK行M列のプロファイル行列（「クロマトグラム行列」ともいう）Ｐの積ＳＰで近似するように、スペクトル行列Ｓ及びプロファイル行列Ｐを求める（ステップ３）。スペクトル行列Ｓは前記第１行列Ｓに相当し、プロファイル行列Ｐは前記第２行列Ｐに相当する。スペクトル行列Ｓ及びプロファイル行列Ｐはそれぞれ、以下のように表される。

【数5】

【数6】

【0042】

ここでスペクトル行列Ｓの各行列要素s_nk（nは1～Nのいずれか、kは1～Kのいずれか）は、試料に含まれるK種の成分の1つ（k番目の成分とする）によるマススペクトルの各m/zにおける強度を示している。また、プロファイル行列Ｐの各行列要素p_kmは、k番目の成分によるクロマトグラムの各保持時間における強度を示している。すなわち、上掲のスペクトル行列Ｓ中において破線で囲んだ行列要素は成分毎のマススペクトルを示し、プロファイル行列Ｐ中において破線で囲んだ行列要素は成分毎のクロマトグラムを示している。スペクトル行列Ｓの各行列要素s_nk及びプロファイル行列Ｐの各行列要素p_kmはいずれも、0又は正の値（非負値）を有する。図３に、３次元データ７１、並びにマススペクトル７２及びクロマトグラム７３のデータの一例を用いて、データ行列Ｘ、並びにスペクトル行列Ｓ及びプロファイル行列Ｐを概念的に図示する。

【0043】

次に、m/z検出部４３は、得られたスペクトル行列Ｓのk番目の列におけるマススペクトルから1又は複数個のピークを検出し、それらピークに対応するm/zを求める、という操作を各列（すなわちkが1～Kの範囲）で行う（ステップ４）。これらのピークに対応するm/zが、目的試料に含まれるk番目の成分に由来するプリカーサイオンのm/zの候補となる。また、保持時間検出部４４は、得られたプロファイル行列Ｐのk番目の行におけるクロマトグラムから1又は複数個のピークを検出し、それらピークに対応する保持時間を求める、という操作を各行（すなわちkが1～Kの範囲）で行う（ステップ５）。これらのピークに対応する保持時間が目的試料に含まれるk番目の成分の保持時間の候補となる。

【0044】

ＭＳ²分析実行条件決定部４５は、これらステップ４で得られたスペクトル行列Ｓのk番目の列におけるm/zの候補と、ステップ５で得られたプロファイル行列Ｐのk番目の行における保持時間の候補より、kが1～Kの範囲内の各々において（すなわち目的試料に含まれる成分毎に）、プリカーサイオンのm/zの候補と保持時間の候補のペアから成るプリカーサイオンリストＬを作成する（ステップ６）。

【0045】

ここで、予め試料を含まないバックグラウンドのＬＣ／ＭＳのデータを取得している場合には、ＭＳ²分析実行条件決定部４５は、当該バックグラウンドのデータに基づいて、プリカーサイオンリストＬからバックグラウンドに由来するm/zの候補と保持時間の候補のペアを除外し、残りの候補のペアを新たなm/zの候補と保持時間の候補のペアとして選択する（ステップ７）操作を行ってもよい。あるいは、この操作の代わりにステップ７において、当該バックグラウンドのデータに基づいて、バックグラウンドの影響を受けていないm/z及び保持時間の範囲内にあるm/zの候補及び保持時間の候補のペアのみを、新たなm/zの候補と保持時間の候補のペアとして選択するようにしてもよい。また、これらステップ７の操作は省略してもよい。

【0046】

ＭＳ²分析実行条件決定部４５は、得られた（ステップ７を実行する場合にはバックグラウンドに由来するm/zの候補と保持時間の候補のペアを除外した後の）プリカーサイオンリストＬに基づいて、前記試料に含まれる成分のプリカーサイオンの選別及び開裂を伴うＭＳ²分析の実行条件（ＭＳ²分析メソッド）を決定する（ステップ８）。ＭＳ²分析では通常、このような分析メソッドは成分毎に既知であることから、そのような既知の分析メソッドを記憶部６１に記憶させておいたうえで、プリカーサイオンリストＬに含まれるm/zや保持時間の候補の情報に基づいて、ＭＳ²分析実行条件決定部４５が記憶部６１から分析メソッドを読み出すようにすればよい。

【0047】

分析メソッドを決定する際に、1回のＭＳ²分析のループタイム（サンプリング間隔）が所定値以下となるように、1つの試料に対するＭＳ²分析の分析メソッドを複数に分割することができる。これにより、十分なサンプリングレートを確保することができ、定量の感度を高くすることができる。

【0048】

また、分析メソッドを決定する際に、分析対象の成分毎に異なるコリジョンエネルギーを設定することができる。これにより、成分毎の適切なコリジョンエネルギーが既知である場合に、成分毎の分析をより適切に行うことができる。あるいは、分析対象の成分毎にそれぞれ複数のコリジョンエネルギーを設定してもよい。これにより、成分毎の適切なコリジョンエネルギーが不明である場合に、複数のコリジョンエネルギーで試行することにより最適な開裂条件を求めることができる。

【0049】

こうして決定された分析メソッドに基づいて、ＬＣ／ＩＴ－ＴＯＦＭＳ１は、ＬＣ／ＭＳ²分析実行部５２がＬＣ／ＩＴ－ＴＯＦＭＳ１の各部を制御することにより、ＬＣ／ＭＳ²分析を実行する（ステップ９）。ＬＣ／ＭＳ²分析は、プリカーサイオンリストＬに含まれる全ての保持時間において実行される。全ての保持時間におけるＬＣ／ＭＳ²分析が完了したとき、ＬＣ／ＩＴ－ＴＯＦＭＳ１の一連の動作が終了する。

【0050】

(3)本実施形態のデータ解析方法（行列分解の操作）の詳細
次に、図４を用いて、本実施形態のデータ解析方法である、行列分解実行部４２において実行される行列分解の操作（ステップ３）の詳細を、行列分解実行部４２内の各部の機能と共に説明する。

【0051】

まず、正則化パラメータ・正則化関数準備部４２１は、複数個（ここではb個とする。bは自然数。）の正則化パラメータの候補λrと1つの正則化関数Ｒ(Ｓ, Ｐ)を用意する（ステップ３１、正則化パラメータ・正則化関数準備工程）。本実施形態では、正則化関数Ｒ(Ｓ, Ｐ)として、行列ＳのL1ノルムと行列ＰのL1ノルムの和であるＲ(Ｓ, Ｐ)=|Ｓ|₁+|Ｐ|₁を用いる。ここで行列のL1ノルムとは、該行列の全ての行列要素の和をいう。正則化パラメータの候補λrは、正の実数を複数個、適宜定める。

【0052】

次に、行列候補決定部４２２は、b個の正則化パラメータの候補λrの各々について、損失関数L(Ｓ, Ｐ) = D(Ｘ|ＳＰ) + λrR(Ｓ, Ｐ) が最小となるように、Ｓrt, Ｐrtを求める。ここで距離関数D(X|SP)=Σ_n,mＤ(Ｘ_nm｜(ＳＰ)_nm)は、該データ行列の行列要素Ｘ_nm及びそれに対応する該ＳＰの行列要素(ＳＰ)_nm毎の距離の総和を示しており、データ行列Ｘと積ＳＰの各行列要素の相違度を示す関数である。本実施形態では要素ごとの距離関数Ｄ(x｜y)として一般化ＫＬダイバージェンスD_KL(x｜y)=x log(x/y)-(x-y)を用いる。

【0053】

さらに、行列候補決定部４２２は、以下の2つの方法のいずれかによって、正則化パラメータの候補λr毎にスペクトル行列候補（第１行列候補）Ｓr及びプロファイル行列候補（第２行列候補）Ｐrと決定する（行列候補決定工程）。第１の方法では、仮候補の組(Ｓrt, Ｐrt)をそのままスペクトル行列候補Ｓr及びプロファイル行列候補Ｐrと決定する（ステップ３３）。

【0054】

第２の方法では、ステップ３３に代えて以下の操作を実行する。まず、Ｓrt, Ｐrtを初期値として、正則化項を含まない第２の損失関数 L₂(Ｓ, Ｐ)= D(Ｘ|ＳＰ) を最小化するように行列Ｓrt2及び行列Ｐrt2を求める（ステップ３３－１）。そして、これら行列Ｓrt2及び行列Ｐrt2をスペクトル行列候補Ｓr及びプロファイル行列候補Ｐrと決定する（ステップ３３－２）。

【0055】

これら2つの方法のいずれかによってスペクトル行列Ｓ及びプロファイル行列Ｐの候補の組(Ｓr, Ｐr)を定めた後、確率分布変換部４２３は、b個の正則化パラメータの候補λrの各々について、距離関数D(X_nm｜(SP)_nm)に対応する確率分布の累積分布関数F_nmを用意する（ステップ３４）。そして、b個の正則化パラメータの候補λrの各々について、データ行列Ｘの行列要素X_nm、及びそれに対応するスペクトル行列候補Ｓrとプロファイル行列候補Ｐrの積ＳrＰrの行列要素(SrPr)_nmの組み合わせ(X_nm, (SrPr)_nm)毎に、それらの行列要素の値を累積分布関数Ｆ_nmに代入することにより、標準一様分布に従うことが期待される y_nm=F_nm(X_nm｜(SrPr)_nm)を求める（ステップ３５）。本実施形態で用いている損失関数である一般化ＫＬダイバージェンスD_KL(Ｘ｜ＳＰ)=Ｘlog(Ｘ/ＳＰ)-(Ｘ-ＳＰ)に対応する累積分布関数Ｆ(Ｘ｜ＳＰ)は、ポアソン分布の累積分布関数であることが知られている。

【0056】

次に、適合度算出部４２４は、b個の正則化パラメータの候補λrの各々について確率分布変換部４２３で求めたy_nm=F_nm(X_nm｜(SrPr)_nm)の経験分布と、標準一様分布の累積分布関数との適合度を算出する（ステップ３６、適合度算出工程）。適合度の算出には、コルモゴロフ－スミルノフ（Kolmogorov-Smirnov：ＫＳ）統計量、クラーメル－フォン ミーゼス（Cramer-von Mises）統計量、又はアンダーソン－ダーリング（Anderson-Darling）統計量と呼ばれる、統計学における既知の適合度計算手法を用いることができる。

【0057】

次に、行列決定部４２５は、正則化パラメータの候補λr毎に算出された適合度の値を対比して、適合度が最大値となる正則化パラメータの候補λrに対応するスペクトル行列候補Ｓr及びプロファイル行列候補Ｐrを、スペクトル行列Ｓ及びプロファイル行列Ｐとして決定する（ステップ３７、行列決定工程）。なお、適合度が最大値となる正則化パラメータの候補λrの代わりに、適合度が所定の閾値以上である正則化パラメータの候補λrのうち、対応するスペクトル行列候補Ｓr及びプロファイル行列候補Ｐrにおける正則化パラメータλrが最大であるものをスペクトル行列Ｓ及びプロファイル行列Ｐとして決定してもよい。以上により、行列分解の操作が完了する。

【0058】

(4) 行列分解の計算の例
次に、図５の左図に示したクロマトグラムと右図に示したマススペクトルのデータを用いて、行列分解の計算を行った例を説明する。ここで、左図では、1回の測定で得られたクロマトグラムを、4種類の成分にそれぞれ起因する合計4つのピーク及びバックグラウンド（BG）という、合計5つに分離して示した。これら5つのクロマトグラムは、データ行列の各行にそれぞれ対応する。従って、この例では、実際の（実験で得られた）データにおけるKの値は5である。右図では、これら5つのクロマトグラムの各々について、対応するマススペクトルを示した。

【0059】

図６に、実際の（実験で得られた）データにつき、m/zが異なるマスクロマトグラムを多数重ねて示す。この図は、データ行列Ｘの情報を示している。すなわち、1つのマスクロマトグラムにおける保持時間毎の強度の値の集合は、データ行列Ｘの1つの行における行列要素の値の集合に相当する。また、或る保持時間における多数のマスクロマトグラムの強度の値の集合は、データ行列Ｘの1つの列における行列要素の値の集合に相当する。

【0060】

このようなデータ行列Ｘに対して、λrが3つの値（1、256、512）の場合についてそれぞれ、本実施形態の行列分解の方法によりスペクトル行列候補Ｓr及びプロファイル行列候補Ｐrを決定し、適合度を求めた。

【0061】

図７に、λr=1の場合について求めた、プロファイル行列候補Ｐrの行列要素の各行をクロマトグラフとして示す（左図）と共に、スペクトル行列候補Ｓrの行列要素の各列をマススペクトルとして示す（右図）。同様に、図８に、λr=512の場合について求めたクロマトグラフとマススペクトルを示し、図９に、λr=256の場合について求めたクロマトグラフとマススペクトルを示す。図７～図９には合わせてＫＳ統計量の計算値を示す。ＫＳ統計量はＫＳ検定により得られる数値であって、この数値が小さいほど適合度が高いことを意味している。

【0062】

λr=1の場合（図７）には、ＫＳ統計量は0.0924である。また、求められたプロファイル行列候補Ｐr及びスペクトル行列候補ＳrにおけるKの値は7となり、実際の値（K=5）よりも大きい。これは、λrが小さ過ぎることにより正則化が不十分であることを意味している。実際、得られたクロマトグラフ及びマススペクトルは、一見して明らかに実際のデータ（図５）と一致していない。

【0063】

λr=512の場合（図８）には、ＫＳ統計量は0.2652である。また、求められたプロファイル行列候補Ｐr及びスペクトル行列候補ＳrにおけるKの値は2となり、実際の値よりも小さい。これは、λrが大き過ぎることにより正則化の効果が強くなり過ぎたことを意味している。得られたクロマトグラフ及びマススペクトルは、一見して明らかに実際のデータ（図５）と一致していない。

【0064】

それに対してλr=256の場合（図９）には、ＫＳ統計量は0.0164であって、3つの候補の中で最小となっている。すなわち、適合度は3つの候補の中で最大である。従って、3つの候補の中からλr=256の場合におけるプロファイル行列候補Ｐr及びスペクトル行列候補Ｓrが、プロファイル行列Ｐ及びスペクトル行列Ｓとして得られる。得られたプロファイル行列Ｐ及びスペクトル行列ＳにおけるKの値は5となり、実際のデータと一致している。また、これらプロファイル行列Ｐ及びスペクトル行列Ｓから求められるクロマトグラフ及びマススペクトルは、実際のデータ（図５）に近いものとなっている。

【0065】

ここでは3つの正則化パラメータの候補λrからそれぞれ得られたプロファイル行列候補Ｐr及びスペクトル行列候補Ｓrによるクロマトグラフ及びマススペクトルを示したが、正則化パラメータの候補λrは3つには限られない。正則化パラメータの候補λrの数を多くするほど、より正確なプロファイル行列Ｐ及びスペクトル行列Ｓを求めることができる。

【0066】

(5) その他
上記実施形態では、データ行列ＸをN行M列、スペクトル行列ＳをN行K列、プロファイル行列ＰをK行M列でそれぞれ規定したが、それらの代わりに、データ行列ＸをM行N列、スペクトル行列ＳをK行N列、プロファイル行列ＰをM行K列でそれぞれ規定してもよい。この場合にはプロファイル行列Ｐが前記第１行列に相当し、スペクトル行列Ｓが前記第２行列に相当する。また、積ＳＰの代わりに積ＰＳ（ここでは、Ｓ及びＰは第１行列及び第２行列の記号ではなく、この例におけるスペクトル行列及びプロファイル行列の記号である）を用いる。

【0067】

クロマトグラフ質量分析装置の構成は上記のＴＯＦＭＳ１には限定されない。例えば、上記実施形態で用いたイオントラップ３０の代わりに、四重極マスフィルタ等のマスフィルタとコリジョンセルを組み合わせたものを用い、上記実施形態で用いたＴＯＦ３１の代わりに直交加速型のＴＯＦを用いたクロマトグラフ質量分析装置に本発明を適用してもよい。また、ＴＯＦＭＳ以外のクロマトグラフ質量分析装置に本発明を適用してもよい。

【0068】

本発明はクロマトグラフ質量分析装置で得られるデータの他に、イメージング質量分析装置やイメージングFT-IR等で得られるデータにも適用することができる。イメージング質量分析装置では、第１の変数及び第２の変数のいずれか一方がm/zであって他方が位置を表す数値に該当する。イメージングFT-IRでは、第１の変数及び第２の変数のいずれか一方が吸収波長であって他方が位置を表す数値に該当する。

【0069】

その他、本発明は上記実施形態には限定されず、種々の変形が可能であることは言うまでもない。

【0070】

［態様］
上述した例示的な実施形態は、以下の態様の具体例であることが当業者により理解される。

【0071】

（第１項）
第１項に係るデータ解析方法は、分析装置から得られた、第１の変数の値が互いに異なるN行のデータと第２の変数の値が互いに異なるM列のデータから成るN行M列のデータ行列Ｘに基づいて、N行K列の第１行列ＳとK行M列の第２行列Ｐの積ＳＰで該データ行列Ｘを近似するように、該第１行列Ｓ及び該第２行列Ｐを求める方法であって、
正則化パラメータの候補λr（rは1からr_maxまでの自然数）を複数個、及び解のスパース性を誘導する正則化関数R(Ｓ, Ｐ)を1個用意する正則化パラメータ・正則化関数準備工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記データ行列Ｘと前記積ＳＰの相違度を示す距離関数D(Ｘ｜ＳＰ) と、正則化パラメータ候補λrと正則化関数R(Ｓ, Ｐ)の積λrR(Ｓ, Ｐ)との和である損失関数 L(Ｓ, Ｐ) =D(Ｘ｜ＳＰ)+λrR(Ｓ, Ｐ)の値が最小となるように求められる行列Ｓrt及び行列Ｐrtを第１行列Ｓの候補Ｓr及び第２行列Ｐの候補Ｐrとして決定するように最適化問題を解く行列候補決定工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記データ行列Ｘの行列要素X_nm及びそれに対応する前記第１行列候補Ｓrと前記第２行列候補Ｐrの積ＳrＰrの行列要素(SrPr)_nmの組み合わせ毎に、前記距離関数D(X_nm｜(SP)_nm)に対応する確率分布P_nmを共通の確率分布 Pcommonへ変数変換する関数F_nmを用いて、共通の確率分布Pcommonへと変数変換された値である変換値y_nm=F_nm(X_nm｜(SrPr)_nm)を求める確率分布変換工程と、
前記複数個の正則化パラメータの候補λrの各々について、前記変換値y_nmと確率分布Pcommonの累積分布関数との適合度を求める適合度算出工程と、
前記複数個の正則化パラメータの候補λrのうち、前記適合度が最大値となるもの、又は前記適合度が所定の閾値以上であって前記第１行列候補Ｓr及び前記第２行列候補Ｐrにおける前記λrの値が最大であるものについて決定した前記第１行列候補Ｓr及び前記第２行列候補Ｐrを、第１行列Ｓ及び第２行列Ｐとして決定する行列決定工程と
を有する。

【0072】

第１項に係るデータ解析方法によれば、試料に含まれる成分の種類の数が未知であっても、該試料に対する測定により得られたデータ行列Ｘに基づいて、適切な、すなわち、該成分の数に近い因子数Kを有する第１行列Ｓ及び第２行列Ｐを求めることができる。

【0073】

（第２項）
第２項に係るデータ解析方法は、第１項に係るデータ解析方法の前記行列候補決定工程において、前記行列Ｓrt及び前記行列Ｐrtを第１行列候補Ｓr及び第２行列候補Ｐrと決定する代わりに、該行列Ｓrt及び該行列Ｐrtを初期値として、正則化項を含まない第２の損失関数 L₂(Ｓ, Ｐ)= D(Ｘ|ＳＰ) を最小化するように行列Ｓrt2及び行列Ｐrt2を求め、該行列Ｓrt2及び該行列Ｐrt2を第１行列候補Ｓr及び第２行列候補Ｐrと決定する。

【0074】

第２項に係るデータ解析方法によれば、損失関数 L(Ｓ, Ｐ) =Ｄ(Ｘ｜ＳＰ)+λrR(Ｓ, Ｐ)の値が最小となるように求めた行列Ｓrt及び行列Ｐrtを初期値として、さらに正則化項を含まない第２の最適化の判定を行うため、正則化による残差の偏りを補正し、より実際のデータに近い第１行列候補Ｓr及び第２行列候補Ｐrを選択することができる。そのため、最終的に求める第１行列Ｓ及び第２行列Ｐの精度を高くすることができる。

【0075】

（第３項）
第３項に係るデータ解析方法は、第１項又は第２項に係るデータ解析方法において、前記正則化関数Ｒ(Ｓ, Ｐ)が、L1ノルム、若しくはL1ノルムとL2ノルムの線形結合、又は、Ｐの各列の総和が1以下になるように解の制約を加えたうえでＳに対してトレースノルム、det|Ｓ^TＳ|、logdet|Ｓ^TＳ+δI|（但し、Iは単位行列、δは正則化関数を制御するハイパーパラメータ）のいずれかを適用したものである。

【0076】

（第４項）
第４項に係るデータ解析方法は、第１項～第３項のいずれか１項に係るデータ解析方法において、前記累積分布関数が、ガウス分布の累積分布関数、ポアソン分布の累積分布関数、指数分布の累積分布関数、Tweedie分布の累積分布関数のいずれかである。

【0077】

（第５項）
第５項に係るデータ解析方法は、第１項～第４項のいずれか１項に係るデータ解析方法において、前記適合度がコルモゴロフ－スミルノフ検定、クラーメル－フォン ミーゼス検定、及びアンダーソン－ダーリング検定のうちのいずれかを用いて求められる。

【0078】

（第６項）
第６項に係るデータ解析方法は、第１項～第４項のいずれか１項に係るデータ解析方法であって、前記行列要素X_nm毎のノイズの分散σ_nm ²が既知である場合において、前記確率分布Pcommonを標準正規分布とし、前記累積分布関数F_nm(X_nm｜(SrPr)_nm)を (X_nm - (SrPr)_nm)/σ_nmとし、前記変換値y_nm が平均0であると想定した場合の不偏分散

【数7】

を用いて表される-|σ_y ² - 1|を前記適合度として用いる。

【0079】

本発明では、第３項で挙げた種々の正則化関数Ｒ(Ｓ, Ｐ)及び第４項で挙げた種々の累積分布関数Ｆ(Ｘ｜ＳＰ)を適宜用いることができる。また、適合度の算出には第５項又は第６項で挙げた、統計学の分野で知られている種々の方法を適宜用いることができる。

【0080】

（第７項）
第７項に係るデータ解析方法は、第１項～第６項のいずれか１項に係るデータ解析方法において、前記分析装置がクロマトグラフ質量分析装置であり、前記第１の変数及び前記第２の変数のいずれか一方がm/zであって他方が保持時間である。

【0081】

（第８項）
第８項に係るデータ解析方法は、第１項～第６項のいずれか１項に係るデータ解析方法において、前記分析装置がイメージング質量分析装置であり、前記第１の変数及び前記第２の変数のいずれか一方がm/zであって他方が位置を表す数値である。

【0082】

（第９項）
第９項に係るデータ解析方法は、第１項～第６項のいずれか１項に係るデータ解析方法において、前記分析装置がイメージングフーリエ変換赤外分光光度計であり、前記第１の変数及び前記第２の変数のいずれか一方が吸収波長であって他方が位置を表す数値である。

【0083】

本発明は、第７項～第９項に記載のように、クロマトグラフ質量分析装置、イメージング質量分析装置及びイメージングフーリエ変換赤外分光光度計で得られるデータの解析に好適に用いることができる。

【符号の説明】

【0084】

１…液体クロマトグラフ／イオントラップ飛行時間型質量分析装置（ＬＣ／ＩＴ－ＴＯＦＭＳ）
１０…液体クロマトグラフ（ＬＣ）部
１１…移動相容器
１２…送液ポンプ
１３…インジェクタ
１４…カラム
２０…質量分析（ＭＳ）部
２１…イオン化室
２２…ＥＳＩノズル
２３…脱溶媒管
２４…第１段中間真空室
２５…第１イオンガイド
２６…スキマー
２７…第２段中間真空室
２８…第２イオンガイド
２９…分析室
３０…イオントラップ
３１…ＴＯＦ
３２…リフレクトロン電極
３３…イオン検出器
３４…Ａ／Ｄ変換器
４０…データ処理部
４１…データ行列作成部
４２…行列分解実行部
４２１…正則化パラメータ・正則化関数準備部
４２２…行列候補決定部
４２３…確率分布変換部
４２４…適合度算出部
４２５…行列決定部
４３…m/z検出部
４４…保持時間検出部
４５…ＭＳ2分析実行条件決定部
５０…分析制御部
５１…ＬＣ／ＭＳ分析実行部
５２…ＬＣ／ＭＳ²分析実行部
６１…記憶部
６２…表示部
６３…操作部
７１…３次元データ
７２…マススペクトル
７３…クロマトグラム

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】