トップ サイトマップ
関連情報
ホーム > 製品&サービス > コンサルティングサービス > カスタムアルゴリズム開発 > カスタムアルゴリズム開発事例:質量分析データ処理向けハイパースペクトルイメージング・ライブラリの開発

カスタムアルゴリズム開発の事例:
質量分析データ処理向けハイパースペクトルイメージング・ライブラリの開発


 質量分析はサンプル物質の化学的構成を解析する際に用いられる方法です。細胞がサンプルの場合はその中の薬物や代謝産物の解析に用いられ、さらに3次元の人工的な固形物でも可能です。これらサンプルは光やイオンに照射されて、散乱されたイオンの飛行時間が検出器で測定されます。これによりイオンの質量電荷比が決定されます。
 質量分析には様々な方法がありますが、最も重要なものは二次イオン質量分析法(SIMS)とマトリックス支援レーザー脱離イオン化法(MALDI)です。これらにおいてはデータはハイパースペクトルとして格納されます:各ピクセルにはデータの全スペクトルが保存され、これらは飛行時間チャネルへ分割され、そのカウント数が各チャネルに記録されます。

質量分析プロセス概念図 SIMS装置

 質量分析イメージングデータは巨大で複雑なものです。また雑音も存在し、特定のイメージング手法を用いるとアーチファクトが含まれる場合があります。こうした事柄が組み合わされて、このデータ解析と可視化を困難なものにしています。
 これに加えて、イメージフュージョン法を用いれば更なる情報を得ることが出来ます。この方法は、複数の異なるスキャンデータを組み合わせて一つのイメージを作成します。下図の例においては、左図のイメージは低空間解像度でスキャンされており、その高解像度ケースが右図に示されています。しかしながら、各図の下に示してある通り、左図は明確な2つのピークを持つ高質量分析解像度を持ち、右図は低解像度になっています。イメージフュージョンの目的は、2つのイメージを組み合わせて高解像質量分析かつ高空間解像度を持つ一つのイメージを生成することです。


低空間解像および高空間解像度のサンプルスキャン

上の低/高空間解像度のサンプルスキャンに対応する飛行時間プロット:それぞれ高/低質量分析解像度を持つ。

 このようなイメージ統合ツールは、衛星画像のような他の分野で既に広く用いられていますが、そのソフトウェアや手法は微小光学解析分野においては有効性に欠ける場合があります。先述した雑音やアーチファクト、データサイズに関する問題が組み合わされて、質量分析におけるイメージフュージョンは計算量が膨大になります。
 「NAG/NPLハイパースペクトルイメージング・ライブラリ」は、高性能データ操作・解析ルーチン群で構成されたツールです。これらは、MALDIあるいはSIMSのrawデータの読み込みからその処理まで、つまりイメージフュージョンの一貫したワークフローを可能にします。これらルーチン群は、マルチコアマシンを有効に活用するように並列化されています。
 このライブラリーのインターフェイスはポリモルフィッククラスで構成され、MATLAB®とC++で利用可能であり、イメージマトリックスや様々な関数がカプセル化されています。ここではMATLAB®を用いたライブラリ利用の実例を示します。

rawデータの読み込み

 SIMSやMALDI装置で生成された様々な出力rawデータをコンフィグレーションファイルに設定するところから始めます。ピーク選択のための"limits"ファイルやピクセルと飛行時間をグループ分けする"bin size"も指定することが出来ます。

configs = '<path to>\input.config'; % The input.config file points to the raw data and allows bins to be specified. Contents of input.config: % coords = <path to>\raw.coords file % limits = <path to>\limits.txt file #for use in peak selection later % tofs = <path to>\raw.tofs file % props = <path to>\properties.txt file % % #Specify bin sizes, here we group 4x4 sets of pixels together and 16 time-of-flight channels: % xbin = 4 % ybin = 4 % tofbin = 16 img = SpecImageFactory(configs);

 オブジェクトimgは、SpecImageクラスのインスタンスです。これは、ピクセル数×チャネル数のサイズのデータ行列を持ちます。データ行列の構築処理はマルチコア上で良好なスケール性を示します。


マルチコアCPUとK80 CPU上でのデータ行列生成のスケール性。最左バーはオリジナルのNPLによるシリアル実行MATLABコードを示す。

 SpecImageを照会すれば、以下の様にデータ行列に関する情報を取得できます。

x = img.getXDim; % Number of pixels in the x direction y = img.getYDim; % Number of pixels in the y direction A = img.getMatrix; % Returns matrix in sparse format m = img.getNumPixels; % Number of pixels in the image matrix n = img.getNumChannels; % Number of channels in the image matrix: thus the data matrix is m x n

メモリーに格納された行列から直接SpecImageを生成することも可能です。

プリプロセス関数

SpecImageクラスは2つのサブクラス、SIMSpecImageおよびMALDISpecImageを持ちます。これらはそれぞれ、SIMSとMALDI装置に特定の機能がカプセル化されています。
 データ解析の前に、その処理により適した形式へデータを変換する必要があります。この前処理はSISMかMALDIのどちらかを利用するかによって異なります。
 SIMSイメージに対しては、Advanced dead-time correctionとスケーリングを行います。Advanced dead-time correction(ADTC)は、SIMS機器でデッドタイムウィンドウのデータを補正する変換です(デッドタイムウィンドウはデータ解析で偽のアーチファクトを導く可能性があります)。ADTCステップに続いて、ADTCの非線形性を補正するために、データ行列に対角行列を前後に乗算してデータ行列をスケールします:

img.advancedDeadTimeCorrection; img.scaleMatrix;

ADTCの詳細は文献:The accuracy and precision of the advanced Poisson dead-time correction and its importance for multivariate analysis of high mass resolution ToF-SIMS data, Bonnie J. Tyler, Surf. Interface Anal. 2014, 46, 581-590、を参照してください。
 サブクラスSIMSpecImageおよびMALDISpecImageの両方で、ピーク選択が可能です。以前にimgを構築した際に、コンフィグレーションファイルにピーク選択制約(limits)を定義しました。この制約を使用して飛行時間チャネルの数を減らし、データ行列サイズを小さくすることができます。ピークが重複している場合、特にビニング後に警告を受け取ることも可能です:

img.peakSelect();

**Warning: 155 pairs of overlapping peak selection bins have been separated

解析の関数

SpecImageクラスの一部として様々なデータ解析関数が可能です。これらの関数を使用すると、データの次元を縮小したり、データの重要な側面を抽出したり、さまざまな方法でデータを分解したりすることができます。ここでは、データ行列のサイズはm×n、mはピクセル数、nは飛行時間チャネル数であると仮定します。


連続ウェーブレット変換

 データ行列内の与えられたピクセル(行)に対して、様々な飛行時間チャネルのカウントをプロットすることができます:


与えられたピクセルに関する飛行時間チャネルのカウント例

連続ウェーブレット変換では、マザーウェーブレットψのスケールかつ伸張されたバージョンで畳み込むことによって、以下の式のf(x)で示されるこの曲線は分解されます:

\begin{eqnarray} C_{s,k}=\int_{-\infty}^{\infty} f(x) \frac{1}{\sqrt{s}} \overline{ψ} \left( \frac{x-k}{s} \right) dx \end{eqnarray}

マザーウェーブレットψは通常、下図のメキシカンハットあるいはハールウェーブレットと言ったコンパクトな台を持つ単純な関数です:

Mexican Hat wavelet Haar wavelet

グラフ的には、以下に見るように、ウェーブレットの選択されたスケールsに対して繰り返し変換して信号f(x)とマッチさせます:


連続ウェーブレット変換

以下の例では、最初に解析で用いるスケールとピクセルを設定しています。ここでは、20ピクセルと4つのバイナリースケールを選択し、zero end extension(これは領域端でのウェーブレットを処理する方法をつたえるものです)と共にメキシカンハットウェーブレットを用います。

for i = 0:19 pixels(i+1) = i; end for i = 1:4 scales(i) = 2^i; end [coeffs] = img.cwt(scales, pixels, 'MexicanHat', 'Zero');

連続ウェーブレット変換はデータ内のピークの解析に用いることが出来ます。この分析の結果、追加のピーク選択制約を指定することができます。これを行うには関数img.setLimitsを用います。


主成分解析

 m×nデータ行列は、n次元空間のm点(1点が1ピクセルに対応)の集合として抽象化されたものと見ることが出来ます。すると各ピクセルは、nを飛行時間チャネル数としたn次元ベクトルとして表現できます。こうしてn軸群上にこれらの点をプロットすることをイメージできます。主成分分析は、点の位置の最大の広がりが第1の軸に沿って生じ、次に第2の軸に沿って最大、などと続くような新しい軸の集合を見つける事が出来ます。
これは2次元では可視化が容易で、以下の赤線が変換後の軸を示します:


2次元データの主成分分析

新しい軸はローディング、ローディングへのデータの射影はスコアと呼ばれます。高次元のデータに対しては、新たな軸の多くは無視されます。それは、これらが小さなデータの広がりしか持たないためです。よってPCAは次元縮小のために用いることが出来て、新たなデータ行列はかなり小さくなります(僅かな数の列のみ残ります)。以下の例では最初の10次元のみを残しています:

[loadings, scores, eigs, tot] = img.pca(10);

ベクトルeigsは、PCA中にデータ行列から構成された共分散行列の固有値を含みます。eigs/totの計算は、新しい次元毎に考慮される分散の割合の尺度を提供します。


t分布型確率的近傍埋め込み法 (t-distributed stochastic neighbor embedding)

 通常t-SNEと略されますが、これは別の次元縮小の方法です。m×nデータ行列を新たなm×k行列(k<<n)へ変換するもので、次元圧縮法(low-dimensional embedding)の1種として知られています。通常はkは2か3を選択され、つまりデータをRGBの値へ変換するのに用いられるのが一般的です。ここで、主成分解析は線形変換(新たな軸は古い軸の線形結合で表される)でしたが、t-SNEは非線形で、より複雑な多項式による変数間の関係が用いられます。
 t-SNEは高次元空間内の近傍の点を一致させようとします。つまりそれらは低次元空間内で接近している事になります。これを行うには、高次元空間および低次元空間における点が互いに接近している確率分布を計算し、確率分布間の差の尺度であるKullback-Liebler発散として知られる量を反復的に最小化します。典型的には、データ行列を約20次元に縮小するためにはt-SNEの前にPCAを実行します。
 t-SNEのパラメータには多くの選択肢があり、計算速度と精度のトレードオフも検討条件となります。これらパラメータについてはソフトウェアに記載されていますが、詳細については文献:Visualizing Data using t-SNE, Laurens van der Maaten and Geoffrey Hinton, Journal of Machine Learning Research, 9(2008) 2579-2605, and Accelerating t-SNE using Tree-Based Algorithms, Laurens van der Maaten, Journal of Machine Learning Research, 15(2014) 1-21、を参照してください。
 下記の例では、データを3次元まで縮小し、前節のPCAの結果から開始します。どの程度の頻度でスクリーンに出力するかを指定するverbosityといった、様々なオプションパラメータが指定されています。

[embedding, kl_divergence, total_iter] = img.tsne(3,'data',scores,'max_iter',100,'verbosity',25);
** Info: Reducing 4096 x 10 matrix to 4096 x 3 using t-SNE. ** Info: Initial Kullback-Liebler divergence = 22.591861. ** Info: t-SNE Iteration 25, Kullback-Liebler divergence = 21.926958. ** Info: t-SNE Iteration 50, Kullback-Liebler divergence = 2.572930. ** Info: t-SNE Iteration 75, Kullback-Liebler divergence = 1.859599. ** Info: t-SNE Iteration 100, Kullback-Liebler divergence = 1.622093. ** Info: Exiting after the maximum 100 iterations were completed. ** Info: Final Kullback-Liebler divergence: 1.622093.

k平均クラスタリング

 k平均クラスタリングは、データ行列に含まれるm個の点をk個のクラスターへ、各クラスター平均との距離が最小になるように分割します。初期クラスター平均は指定も出来ます。また、標準的な幾何距離かベクトル間の距離算出にコサイン類似度のどちらかを使うことも可能です。以下の例では、SpecImageコンストラクタを用いてPCAスコアから新たなイメージを生成して、k=4としてデフォルトの幾何距離とランダム初期クラスター平均を用いたk平均クラスタリングを実行します。ただしSpecImageコンストラクタを用いるにはスコア行列の転置(')が必要です。

img2 = SpecImage(scores'); [cmeans, inc, nic, wcs] = img2.kmeans(4);

非負値行列因子分解

 非負値行列因子分解は、その積が入力行列Aを近似する2つの行列WとHを算出する手法です。この時、Aは無論、WとHの行列要素は非負値で、その次元はそれぞれm×k、k×nです。ここでkはユーザが指定する小さな次元数です:


非負値行列因子分解は、一つの行列を短い次元と長い次元数を持つ縦棒と横棒のような2つの行列に分解する

 行列Wは基底行列(features matrix)、Hは係数行列とも呼ばれます。以下の例ではk=50として非負値行列因子分解を実行します。

[W,H] = img.nmf(50);

結果

 以下のプロットは、解析ルーチンからの出力を示しています。最初の2つは主成分解析、次にt-SNEを使用した3次元への縮小、4つのクラスターによるk平均クラスタリング、および非負値行列因子分解の出力です。


NAG/NPLデータ解析ルーチンのプロット

より詳しい情報について

 NAG/NPLハイパースペクトルイメージング・ライブラリは、英国国立物理学研究所:National Physical Laboratory (NPL),National Centre of Excellence in Mass Spectrometry Imaging (NiCE-MSI)との2年間にわたる共同作業の成果です。このプロジェクトはInnovate UKの資金で実施されました。
 このライブラリーの計算カーネルはC言語で記述され、一部はNAG LibraryのFortranコードが用いられており、並列処理にはOpenMPが使用されています。最上位レベルはC++およびMATLABのクラスで記述されています。
 このライブラリーには包括的なドキュメント、C++とMATLABによる例題が備わり、MATLABのヘルプシステムに全て組み込まれています。C++に関する文書はonlineで閲覧できます。
 また、64-bit Windows向けのビルドが可能で、24コアで良好なスケール性を示しています。これらのソフトウェアはNPL's internal image fusion softwareにも組み込まれました。

 さらに詳細をお知りになりたい方はコンサルティングサービス・ご相談フォームにご記入の上お問い合わせください。



Results matter. Trust NAG.

Privacy Policy | Trademarks