JSBi Bioinformatics Review
Online ISSN : 2435-7022
総説
行列・テンソル分解によるヘテロバイオデータ統合解析の数理
―第4回 質的データ、距離、グラフ―
露崎 弘毅
著者情報
ジャーナル オープンアクセス HTML

2022 年 3 巻 2 号 p. 33-46

詳細
Abstract

生命科学分野で取得されるデータ集合は、雑多(ヘテロ)な構造になり、ヘテロなデータ構造を扱える理論的な枠組みがもとめられている。本連載では、汎用的なヘテロバイオデータの解析手法である行列・テンソル分解を紹介していく。第4回では、これまで扱ってこなかった質的データ、距離、グラフといった特殊なデータに対して適用できる行列・テンソル分解を紹介する。

質的データとは

本連載[1, 2, 3]でこれまで扱ってきたデータは、マイクロアレイの蛍光強度など、最初から数値として計測されるものを対象としていた。このようなデータは量的データという。量的データにはさらに、原点が存在し、間隔や比率に意味がある比率尺度(ratio scale)と、原点は存在しないが、間隔(または値の差)には意味がある間隔尺度(interval scale)とに分類される[4, 5]。例えば、遺伝子発現量やタンパク質発現量など、物質の量を計測したデータは、絶対的な原点(物質量が0)が存在し、200個の物質は100個の物質よりも量が2倍あるという比率が議論できるため比率尺度となる。一方、時間や摂氏温度などは、便宜上0となる点は設定できるものの、そのような0に本質的な原点の意味は無い。そのため、温度が5度から10度になったからといって、実質的に2倍の何かが増えたわけではないため、間隔尺度となる。

一方、次世代シーケンサーや一塩基多型(Single Nucleotide Polymorphism;SNP)を計測するマイクロアレイから出力された塩基配列情報や、遺伝子とアノテーション情報の対応表など、数値では表現されないデータもある。このようなデータは質的データという。質的データにはさらに、データ間に大小関係が存在する順序尺度(ordinal scale)と、大小関係が存在しない名義尺度(nominal scale)とに分類される[1, 2]。例えば、疾患の重症度(ステージ)は、数値データではないが、ステージが進むほど疾患の症状が悪化するため順序が存在する、すなわち順序尺度である。名義尺度とは、DNAを構成する塩基(A/T/G/C)や、異なる実験条件(例:試薬投与群/対照群)など、単に区別があるだけで、値同士に順序は存在しないデータである。

質的データは量的データよりも扱いが難しい。というのは、数値ではないため、量的データでは当たり前にできる1+2=3のような数値同士の計算が、質的データになるとできないためである。例えば、名義尺度同士のA+Tのような演算は自明ではない。適当にA→1、T→2、G→3、C→4と数値に対応させても、量的データとして扱っていいことにはならない(アンケート情報など、順序尺度では広く利用されるやり方、Likert尺度という)。そのため、これまで本連載で紹介してきたような量的データを想定している行列・テンソル分解手法を、直接質的データに適用することはできず、何らかの工夫が必要となる。

質的データを、量的データに変換するやり方として、一つはなんらかの方法でカウントデータに変換して扱う方法がある(図1)。例えば、SNPアレイのデータ解析では、リファレンス配列と比較して、ある個人のSNPの両側アリルともリファレンス側と同じであった場合(例:AA)、片側アリルだけ変異が入った場合(例:AC)、両側アリルとも変異が入った場合(例:CC)とで、各々0、1、2という数値に変換する。この操作はCの数、またはリファレンスからの変異の蓄積度合い、逸脱度を定量化しているため、{0, 1, 2}のどれかで値が埋められた個人×SNPのカウント行列に対して、主成分分析(Principal Component Analysis;PCA)が適用される[6, 7, 8]。また、RNA-Seqのデータ解析では、既知の遺伝子領域にリードが幾つマッピングされたのかを数えることで、遺伝子ごとのカウント値にしてから、各種解析手法が適用される[9]。ただし、上記のようなやり方は、AからCへの変異など、変化の方向を事前に定める必要があったり、比較するための事前情報(例:リファレンス配列、遺伝子領域)が必要であることから、どのようなデータでも適用できる方法ではない。

図1:質的データの量的データへの変換

SNPの場合、リファレンス配列との比較によるカウント行列化や、取りうる両側アリルのパターンを列挙したダミー変数行列が考えられる。

より汎用的な方法としては、ダミー変数の利用が挙げられる(図1)。ダミー変数は、機械学習の分野では、one-hotベクトルとも呼ばれるもので、ある質的データに対して、とりうる値を全て列挙し、そのどれかに該当した場合は1、しない場合は0と記述することで、二値ベクトルにする方法である。例えば、上記のSNPの話しであれば、あるSNPサイトの変異の仕方は、アリルごとに{A, T, G, C}の4パターンあるため、両アリル分で4×4=16パターンある。あるSNPがACであった場合には、この長さ16のベクトルのうち、AC部分の1要素だけが1、それ以外の15要素は0が埋められる(図1)。このダミー変数ベクトルを複数SNP分まとめることで、ベクトルのオーダーが1つ上がり、ダミー変数行列となり、多変量解析の手法が適用できるようになる。

なお、ダミー変数行列を中心とした多変量解析は、数量化理論[4]として体系づけられており、ダミー変数行列における重回帰分析は数量化Ⅰ類、Fisherの判別分析は数量化Ⅱ類、正準相関分析(Canonical Correlation Analysis;CCA、第2回[2]を参照)は数量化Ⅲ類、後述する多次元尺度構成法(Multi-Dimensional Scaling;MDS)は数量化Ⅳ類と呼ばれる。

分割表の仮説検定

ここでは、2種類のダミー変数行列Zrk×n)とZck×m)があるとする(図2)。この時、これらダミー変数行列同士の行列積(n×m)   

(1) O = Z r T Z c
は、分割表(またはクロス集計表)といい、n通りの行変数とm通りの列変数のとりうる全ての組み合わせごとに、カウント値を集計した結果となる。ここで、分割表の行ごとに和をとったものを、行の周辺度数Oi.という(図2)。なお、分割表の全要素の和は   
(2) N = i = 1 n j = 1 m O i j
であることから、Oi.の総和は   
(3) i = 1 n O i . = N
となる。Oi.をさらにNで割ったものは、行の周辺確率Pi.という(図2)。

  

(4) i = 1 n O i . N = i = 1 n P i . = 1

図2:分割表からもとめられる各種確率

行変数・列変数のダミー変数行列の行列積により得られる分割表を、さらに同時確率行列や、行・列の周辺確率へ変換している。

同様に、列ごとに和をとったものを、列の周辺度数O.j、さらにO.jNで割ったものを、列の周辺確率P.jという(図2)。

  

(5) j = 1 m O . j = N

  

(6) j = 1 m O . j N = j = 1 m P . j = 1

ここで、行変数と列変数とが互いに独立であると考えた場合、ij列目の要素の確率は、Pi.P.jの積として表される。

  

(7) P ( i , j ) = P i . P . j

そのため、分割表のij列目の値の期待値は、   

(8) E i j = N P i . P . j
となる。なお実データの解析においては、行変数と列変数が完全に独立であることはほとんど無いので、観測値と期待値がどの程度ずれているのかを定量的に評価するために、以下のχ2統計量が利用される。

  

(9) χ 2 = i = 1 n j = 1 m ( O i j E i j ) 2 E i j = i = 1 n j = 1 m ( O i j N P i . P . j ) 2 N P i . P . j

この値が大きいほど、観測値と期待値にずれがある、すなわち行変数と列変数がより独立ではない状態にあるということになるが、そのずれが偶然とするならばどの程度起きづらいことなのかを調べるためにχ2検定が利用される。χ2検定では、χ2統計量が自由度(n-1)×(m-1)のχ2分布に従うとして、実データから求められた統計量よりも裾野にある部分の面積(p値)を計算し、解析者が事前に設定した有意水準(多くの場合、0.05や0.01)よりも、p値が小さいか否かを調べる。行変数と列変数の独立性検定の仕方としては、他にも尤度比検定(G検定)、Fisherの正確確率検定などが挙げられる[4]。

分割表の分解

上記の仮説検定を分割表に適用することにより、ある分割表の行変数と、列変数が独立ではないか(またはどちらとも言えないか)を調べることはできる。ただし、仮説検定が示すことは、あくまで分割表全体の傾向であり、より具体的に、何行目と何列目の関連性が、独立性を妨げているのかにまでは踏み込めない。対応分析(Correspondence Analysis;CA)は、分割表に行列分解を適用することで、データからそのような情報を読み取るための手法である[4, 10, 11, 12]。

ここでは再び上記のn×mの分割表を考える。CAでは以下のPearson残差(または標準残差)という尺度を利用する。

  

(10) χ 2 N = i = 1 n j = 1 m P i j P i . P . j P i . P . j

ただし、PijOijNは同時確率行列である(図2)。行方向の周辺確率をr(=P1n、1nは長さnで全ての要素が1のベクトル)、列方向の周辺確率をc(=PT1m、1mは長さmで全ての要素が1のベクトル)として、rを対角要素に並べた対角行列をDrcを対角要素に並べた対角行列をDcとすると、分割表の個々の要素のPearson残差は以下のように行列表記でまとめて表すことができる。

  

(11) S = D r 1 2 ( P r c T ) D c 1 2

CAは、この行列に対して特異値分解(Singular Value Decomposition;SVD、第1回[1]を参照)を適用することで、Sを以下のように分解する。

  

(12) S = U Λ V T

式(10)よりPearson残差の自乗はχ2Nであるため、   

(13) χ 2 N = trace ( S T S ) = trace ( V Λ U T U Λ V T ) = trace ( V Λ 2 V T ) = trace ( Λ 2 V T V ) = trace ( Λ 2 )
となり、分割表の行変数・列変数の独立性(χ2)は、行列Sの特異値の自乗和/STSの固有値の和(Λ2)で決定されることがわかる。なおSVDの圧縮次元数kは、最大でmin(n, m)-1まで設定可能ではあるが、それより小さい次元で打ち切った場合、CAχ2統計量の最小近似を提供するようなUV、言い換えれば、互いに非独立で関連性の高い、行変数のスコアと列変数のスコアを同時に探索することに相当する。具体的に、どの行変数とどの列変数で関連性が高いのかは、U, Vを同じ空間に射影するバイプロットで確認できる。

なお、実際のCAでは、重み付きSVD(weighted SVD13])と呼ばれる特殊なSVDが利用されており、通常のSVDとは拘束条件に違いがある。例えば、通常のSVDの拘束条件は、UVの列直交性、すなわち以下のものとなるが、   

(14) U T U = V T V = I k
CAでは、   
(15) U ~ T D r U ~ = V ~ T D c V ~ = I k
または、   
(16) U ^ T D r U ^ = V ^ T D c V ^ = Λ
という拘束条件を設定し、U, Vではなく、上記の U ~ , V ~ (標準座標;Standard coordinates)、または U ^ , V ^ (主成分座標;Principal coordinates)を行・列のスコアとして求める。これらのスコアは、以下のように通常のSVDU, Vから求めることができる。

  

(17) U ~ = D r 1 2 U

  

(18) U ^ = D r 1 2 U Λ = U ~ Λ

  

(19) V ~ = D c 1 2 V

  

(20) V ^ = D c 1 2 V Λ = V ~ Λ

Pearson残差は、1細胞RNA-Seqデータの正規化手法として、近年注目されており、細胞ごとのリード数のばらつきにロバストであるという報告がある[14, 15, 16, 17]。CAは、生命科学分野の中でも、群集生態学でいち早く取り入れられてきた[18]。データに並び順が存在し、その並び順に従いデータが連続的に変化する効果が支配的である場合に、1番目と2番目のスコアでプロットすると、U字型の軌道がアーティファクトとして出てしまう「馬蹄効果」が、PCAと比べると比較的起きづらいという理由から、好まれていたようである。バイオインフォマティクス分野でも、生物種ごとのコドン使用頻度[19, 20]、タンパク質のアミノ酸組成[21, 22]、SNP23]など、様々な質的データで利用されている。PCAと比べて、若干解析結果が向上する場合があるとのことで、量的データであっても、あえて質的データに変換してからCAを適用した例もある[24, 25, 26]。また、複数バッチの1細胞RNA-Seqデータをまとめて次元圧縮する際に、バッチごとの違いがデータに大きく表れてしまう「バッチ効果」や、他よりも大きく値が異なる「外れ値」細胞の影響に対して、PCAよりも比較的ロバストであるという報告もある[27]。式(1)を見てわかるように、CAは2つのダミー変数行列の同時分解ととらえることもできるため、第2回[2]で紹介したCCAのように、サンプル×遺伝子発現量行列と、遺伝子×遺伝子オントロジー行列の統合解析として利用された例がある[28]。

多元分割表の分解

行変数と列変数の関係性を扱う上記の二元分割表は、3つ以上の変数を扱えるように拡張でき、その場合の分割表は多元分割表という。多元分割表に対する古典的な仮説検定手法としては、K個の2×2分割表(K≥2)において、全ての分割表のオッズ比が等しいことを帰無仮説とするCochran-Mantel-Haenszel検定(CMH検定)[4]や、行変数と列変数と3つ目の変数(奥行き変数)の中から2つの変数を選び、どの変数間に交互作用があるのかを調べる対数線形モデル[4]といった手法がある。しかしながら、これら仮説検定を利用しても、多元分割表全体の傾向はわかるが、具体的に多元分割表内のどの位置の値が独立性を妨げているのかまではわからないのは、二元分割表の時と同様である。

CAは、多元分割表にも拡張できる。一つのやり方としては、多重対応分析(Multiple Correspondence Analysis;MCA29])がある。MCAでは、属性が異なるK個の二元分割表を連結させて、1つの二元分割表にしてからCAを適用する。しかしながら、連結する際に、属性の違いはわからなくなってしまうデメリットがある。二元分割表に変形せずに、多元分割表をそのまま扱う手法としては、多元分割表がデータ構造としては高次のテンソルであることに着目し、CP分解やTucker分解など、一般的なテンソル分解手法を適用したThree-way CAという手法がある[30, 31, 32, 33]。ただし、生命科学分野での利用例は筆者の知る限りまだない。

距離とは

2つのデータが似ている度合いを表現する尺度として、類似度と非類似度がある[5, 34]。類似度は、値が大きいほど2つのデータが似ていることを意味し、非類似度は反対に、値が小さいほど2つのデータが似ていることを意味する。代表的な類似度としては、本連載[1, 2, 3]で何度も扱ってきた内積、共分散、相関係数などが挙げられる。非類似度の中でも特に、データa, b, cに関して以下の距離の公理を満たす尺度dのことを距離(または計量)という。

  

(21) d ( a , b ) 0 ( X 1 ) d ( a , b ) = 0 X = Y ( X 2 ) d ( a , b ) = d ( a , b ) ( X 3 ) d ( a , b ) + d ( b , c ) d ( a , c ) ( X 4 )

条件(X1)は距離が負にならないこと、条件(X2)は自分自身との距離は0になること、条件(X3)は、距離が対称であることを意味する。条件(X4)だけはやや難しいが、三角不等式と呼ばれるもので、3点a, b, cが描く三角形において、任意の2辺の長さの和が残りの1辺の長さ以上になることを意味する。公理は数学の理論を構築する上での前提であり、証明するものではないが、どれも自然な前提だと思われる。

距離の公理を満たすものとして、広く利用されるユークリッド距離の他にも、ハミング距離、マンハッタン距離、チェビシェフ距離、ミンコフスキー距離、マハラノビス距離、レーベンシュタイン距離など、実に多くの距離尺度が提案されている[5, 34]。距離の公理を満たさない非類似度もある。例えば、確率分布間の非類似度を示すKullback-Leibler Divergence(KL距離)は条件(X3)を満たさず、正式には距離ではないため、疑似距離と呼ばれる。

今、2つのベクトルa, b間のユークリッド距離の2乗を考えると、   

(22) a b 2 = a 2 + b 2 2 a T b
となり、以下のような形に直すことができる[35, 36, 37, 38]。

  

(23) a T b = 1 2 ( a b 2 a 2 b 2 )

この式は、内積aTbは、abのベクトルの長さ(L2ノルム)や、a, b間のユークリッド距離から求められることを意味する。このことは、すぐ後で紹介する古典的MDSを考える上で非常に重要な性質である。

距離行列の分解

ここでは、n個のデータ間の総当たりの距離が格納された、n×nの距離行列Dを考える。n×nの類似度行列の場合は、何らかの方法で距離の公理を満たすように変換を行う。例えば、n×nのピアソン相関係数行列(-1≤C≤1)を距離に変換するために、D=1-absC)や、D=1-CC(∘は行列間のアダマール積)といった変換がよくなされる。他にも、類似度⇄距離・非類似度の変換の仕方は、数多く存在するので、R言語のsmacofパッケージのsim2diss関数などを参照して欲しい[39]。

古典的MDS、または主座標分析(Principal Coordinate Analysis;PCoA)とも呼ばれる手法では、距離行列Dが与えられている時に、データの座標Yn×k)を復元する手法である[35, 36, 37, 38](図3)。式(22)のユークリッド距離と内積の関係性を、n個のデータに拡張して行列表記にすると、距離行列の2乗は、   

(24) D 2 = g 1 n T + 1 n g T 2 Y Y T
となる。ただしgは、内積行列(グラム行列)YYTの対角要素を含んだ長さnのベクトル(gdiagYYT))、1nは長さnで全ての要素が1のベクトルである。そして、式(23)と同様、左辺をグラム行列に直すと、   
(25) Y Y T = 1 2 H D 2 H
となる。ただし H = I 1 n 1 n 1 n T In×nの単位行列である。n×nの行列に作用することで、その行列の平均値を0にすることから、Hは中心化行列(Centering Matrix)と呼ばれる。Hを左右からかける操作は、D2の行方向と列方向とで、ともに平均を0にするため、Double Centering、またはYoung-Householder変換[35, 36, 37, 38]と呼ばれる。なお、式(25)の左辺のYは、この段階ではまだわかっていないが、グラム行列GYYTを固有値分解すると、   
(26) G = U Λ U T
となり、固有ベクトルを列ごとに格納したn×kの行列Uと、対角要素ごとに固有値を格納したk×kの対角行列Λに分解できることから、Y1/2としてもとめられる。なお、圧縮次元数kは、最大でn-1まで設定可能ではあるが、それより小さい次元で打ち切った場合、Yは元の次元でのデータ間の距離関係をなるべく保持するような低次元座標としてもとめられる。

図3:MDSのテンソル拡張

1つの距離行列に適用するMDSに対し、INDSCALでは複数の距離行列をまとめた3階テンソルに対してテンソル分解を適用する。

古典的MDSは、SNPなどゲノム配列の特徴による集団構造の推定[40, 41, 42, 43]、マルチプルアライメントの結果を利用した配列の分布推定[44]、遺伝子発現量による臓器・サンプルの分布推定[45, 46, 47]、メタゲノムデータの16Sプロファイルの分布推定[48]など、バイオインフォマティクス研究での応用例は多い。近年だと、Hi-Cで得られたコンタクトマップから染色体の三次元構造を推定する問題にも応用されている[49, 50, 51]。古典的MDSは、データ行列を一度n×nの距離行列に変換してから適用される場合も多いが、データの性質や計測の仕方の都合上、個々のデータのベクトル表現は不明だがデータ間の類似度・非類似度だけはわかる状況、すなわち「最初からn×nの行列しか与えられていない場合」(例:Hi-C)にも適用可能なのも魅力である。後者の場合では、PCAや独立性分分析(Independent Component Analysis;ICA)、非負値行列分解(Non-negative Matrix Factorization;NMF)といったデータ行列に対する行列分解手法はそのままでは適用できないことから、古典的MDSが最初に利用される手法となりえる。

古典的MDSから派生して、データ間の距離ではなく、データ間の非類似度関係の順序だけ一致していれば良いと仮定することで、非線形な次元圧縮を行う非計量MDSや、複素数にまで拡張することで、非対称な行列でもMDSで扱えるようにした非対称MDSといった手法も提案されている[35, 36, 37, 38]。他にもMDSは非線形な次元圧縮法を扱う多様体学習[52]という分野に発展しており、古典的MDSの距離行列Dに重み行列WをかけたSammon Mapping、k近傍のデータとの距離関係(測地線距離)だけを利用してMDSを行うISOMAPやLLE、高次元・低次元空間での距離行列の計算に非線形な射影を行うt-SNEやUMAP[53]といった手法が次々と提案されている。特に、t-SNEとUMAPは近年1細胞ゲノミクス分野で広く利用されている。また、非計量MDSに関しても、ニューラルネットワークなどの非線形な次元圧縮法に拡張したTriplet Embedding(またはOrdinal Embedding)という分野に発展しており[54]、t-SNEやUMAPと類似したt-STE[55]やTriMAP[56]といった手法が登場している。

複数の距離行列の同時分解

古典的MDSのテンソル拡張として、Individual Differences Scaling(INDSCAL35, 36, 37, 38])がある(図3)。INDSCALでは、l個のn×n距離行列D1, D2, ..., Dlから、l個のn×nグラム行列G1, G2, ..., Glを計算し、全距離行列間で共通する行列Un×k)と、距離行列ごとに得られる対角行列Λ1, Λ2, ..., Λlk×k)を以下のように推定する。

  

(27) G o = U Λ o U T ( o = 1 l )

1個体ごとに1距離行列が得られる状況では、距離行列ごとの重みの違いを仮定しない場合、例えば、複数のグラム行列の平均行列GaverageGaverageUΛUTとして分解した場合と比べてみると、個体ごとに異なるΛoによって、Uの列ベクトルのスケールを個体ごとに調整できるため、INDSCALは個体差MDSとも呼ばれる。l個のn×nグラム行列は、奥行き方向に格納していくことで、n×n×lの三階テンソルになる。そして、INDSCALは、第3回[3]で紹介したテンソル分解の一つCP分解の特殊な場合とみなせるため、UΛoはCP分解の解法である、Alternating Least Squares(ALS)により逐次的に推定される[35, 36, 37, 38]。一方、ALSによる解法は、非対称距離行列に適用できない、重みが負の値をとりうる、低次元座標が互いに似通ってしまうdegeneracy(縮退)という現象が起きる、といったことが問題視されており、SVDを利用して互いに直交な低次元座標を求めるOrthogonal INDSCAL57]という手法も提案されている。

バイオインフォマティクス分野におけるINDSCALの適用事例はまだそれほど多くないが、メタボローム解析で幾つか利用されており、サンプルの違いを考慮して、代謝物の分布を可視化できる有効性が示されている[58, 59, 60, 61]。また、ケモインフォマティクス分野で、化合物の分布の可視化にも利用されたことがある[62]。

グラフとは

関数やデータ点の分布などの情報を視覚的に表現した図(プロット、フィギュア)のことをグラフと言う場合もあるためややこしいが、ここではグラフ理論など、ネットワークの意味でのグラフを扱うこととする。グラフ[63, 64, 65, 66]とは、ノード(頂点)Vと、エッジ(辺)Eの集合で、以下のように定義される。

  

(28) G = ( V , E )

例えば、図4の例では、5つのノードと5つのエッジで構成されたグラフを示しており、このグラフは、   

(29) V = { v 1 , v 2 , v 3 , v 4 , v 5 } E = { e 1 , e 2 , , e 5 } = { { v 1 , v 2 } , { v 2 , v 3 } , { v 2 , v 4 } , { v 3 , v 4 } , { v 3 , v 5 } }
と定義される。以降は、Vの要素数は|V|、Eの要素数は|E|とする。なお、エッジに向きがないものは無向グラフ、向きがあるものは有向グラフと呼ばれるが、本項では無向グラフのみを扱うものとする。

図4:グラフの行列表現

グラフを行列分解の枠組みに載せるために、グラフを行列として扱う。

以降で、グラフを行列分解の枠組みに載せるための準備として、グラフを行列として扱うことを考える。そのための方法として代表的なものを2つ紹介する。一つは、隣接行列(Adjacency Matrix)と呼ばれるもので、ノードviとノードvj間にエッジがあれば1、無ければ0とする|V|×|V|の行列Aとして、以下のように定義する(図4)。

  

(30) A ( i , j ) = { 1 ( { v i , v j } E ) 0 ( otherwise )

ただし、自分自身とのエッジ(自己ループ)は考慮しないため、対角要素は全て0とする。もう一つは、接続行列(Incidence MatrixまたはConnectivity Matrix)と呼ばれるもので、ノードviがエッジejと接続していれば1、していなければ0とする|V|×|E|の行列Hgとして、以下のように定義する(図4)。

  

(31) H g ( i , j ) = { 1 ( v i e j ) 0 ( otherwise )

なお、隣接行列Aと接続行列Hgには、   

(32) A = H g H g T D v
という関係性がある。ただしDvは、i番目の対角要素にノードviに接続するエッジの数(ノードviの次数(degree))、非対角要素には0を持つ次数行列であり、Aの対角要素を0にするためにHgHgTから引いている。また、エッジに重みがあるグラフにおける隣接行列をAweight(対角要素は0、それ以外の要素は非負値)とすると、対応する次数行列Dv, weightと、i番目の対角要素にエッジeiの重みを持つ重み対角行列Wとの間には、以下の関係が成り立つ。

  

(33) A weight = H g W H g T D v , weight

グラフラプラシアンの分解

あるグラフに属するノード集合を2つの部分グラフに分割することを考える。ただし、異なる部分グラフ同士のノード間に引かれたエッジは切断するものとする。エッジの切断の仕方は幾通りも考えられるが、ここではできるだけ少ないエッジを切断することで、2つの部分グラフに分割したい。このような問題は、グラフ理論における最小カット問題となる[63, 64, 65, 66]。

あるノードが2つの部分グラフ(ここでは部分グラフ1、2とする)のどちらに属するのかを示す、長さ|V|のベクトルsを考える。ただしsの各要素は部分グラフ1に属する場合は1、部分グラフ2に属する場合は-1とする。s同士の外積ssTは、|V|×|V|の行列となり、この行列の(i, j)成分は、ノードviとノードvjが同じ部分グラフ同士であれば1、異なる部分グラフ同士であれば-1を示す。そのため以下の行列   

(34) Z = 1 2 ( 1 s s T )
i, j成分は、ノードviとノードvjが同じ部分グラフ同士であれば0、異なる部分グラフ同士であれば1となる。この行列Zと隣接行列Aのアダマール積(∘)の全要素の和の1/2は、以下のようになる。

  

(35) 1 2 i n j m ( Z A ) i , j = 1 4 i n j m ( ( 1 s s T ) A ) i , j = 1 4 i n j m ( A s s T A ) i , j = 1 4 ( s T D v s s T A s ) = 1 4 s T L s

これは異なる部分グラフ同士のノード間に引かれたエッジの本数を意味する。ただし、ここでは1|V|を長さ|V|の全ての要素が1のベクトル、Dsを対角要素にsを格納した対角行列として、 i , j A = i , j D v = 1 | V | T D v 1 | V | = s T D v s i , j s s T A = i , j D s A D s = s T A s と式変形している。最終的に得られた行列   

(36) L = D v A
はグラフラプラシアンと言う。ここでsが取りうる値は、2|V|通りあることから、|V|が大きくなるに従い、組み合わせ最適化の観点からは、sT Lsを最小化するsを厳密に求めることは難しくなる。

一方で、   

(37) s T L s = λ
を以下のように式変形すると、   
(38) L s = λ s
となることから、sT Lsを最小化させるsは、グラフラプラシアンLの最小固有値λに対応する固有ベクトルsで代用できることが期待される。実際には、Lの固有値分解によって厳密に{1, -1}を値に持つベクトルsが得られるわけではないが、類似したパターンを固有ベクトルとして得ることは可能である。このように整数条件を緩和し、変数が連続値をとることを許容した定式化を連続緩和という[66]。

なお、Lの固有値分解で得られる最小の固有値は0であり、その固有値に対応する固有ベクトルでは、全ての要素が定数となることがわかっている。なぜなら、Lの定義からLの行ごとの和は全て0になるが、和をとることは全ての要素が1のベクトル(1ベクトル)をかけることと同じであり、上記の定数ベクトルは、1ベクトルの定数倍だからである。この自明の固有ベクトルは、しかしながら、全てのノードが同じクラスタに属することを意味し、実用上は有用では無いため、2番目に固有値が小さい固有ベクトル(Fiedlerベクトルという[64])から順に注目される。3番目に固有値が小さい固有ベクトルは、2番目の固有ベクトルで分割された後のグラフをさらに分割する。グラフをk個の部分グラフに分割したい場合、小さい順に2番目からk番目までの固有値・固有ベクトルを利用する。得られた固有ベクトルは、グラフ構造を反映させた座標となることから、これを次元圧縮とした場合は、ラプラシアン固有マップ、さらにそれら固有ベクトルに対してk-means法などクラスタリング手法を適用したものは、スぺクトラルクラスタリングという[63, 64, 65, 66, 67]。

正規化の仕方によって、グラフラプラシアンには、様々なバリエーションが存在する。例えば、上記の正規化されていないグラフラプラシアンの定義では、次数が大きいノード同士がまとまって1つの固有ベクトル内で値が大きくなりがちで、その結果スペクトラルクラスタリングとしては、ノードを多く含んだ1つの大きなクラスターと、その他小さいクラスターに分割されやすい。(対称)正規化グラフラプラシアンでは、このバイアスを補正するために、ノード次数によって、以下のようにLを正規化する[63, 64, 65, 66, 67]。

  

(39) L sym = D v 1 2 L D v 1 2 = I D v 1 2 A D v 1 2

著者の知る限りでは、素のLよりもLsymを利用した論文の方が多い印象がある。なお、Lを利用したスペクトラルクラスタリングはレイシオカット(RatioCut)、Lsymを利用したものは、ノーマライズドカット(NCut)として区別される[63, 64, 65, 66, 67]。

また、ランダムウォーク正規化グラフラプラシアンは、   

(40) L rw = D v 1 L = I D v 1 A
として定義されるもので、ガウスカーネルで計算されたグラム行列を(重み付き)隣接グラフとみなして計算したLrwの固有ベクトルは、拡散マップ(Diffusion Map)という次元圧縮法になる[63, 64, 65, 66, 67]。拡散マップは、細胞の分化など、連続的な現象を表現しやすいということで、1細胞RNA-Seq分野で広く利用されている[68, 69, 70, 71]。

他にもラプラシアン固有マップはGraph Embedding(またはNetwork Embedding)[72, 73]という分野に発展しており、深層学習モデルを利用した各種非線型な手法が次々と提案されている。一例として、グラフ上でのランダムウォークで生成した系列データセットに対してword2vecを適用するDeepWalkは、生命医科学RDFデータを要約するのに用いられたことがある[74]。

ハイパーグラフラプラシアン/ラプラシアンテンソルの分解

グラフを拡張した概念として、ハイパーグラフがある。グラフでは、1つのエッジは2つのノードを繋ぐものであったが、ハイパーグラフでは、1つのエッジがあらゆる数のノードを内包するものであり、そのようなエッジはハイパーエッジと呼ばれる[75, 76, 77, 78](図5)。

図5:ラプラシアン固有マップのハイパーグラフ、テンソル拡張

ラプラシアン固有マップは1つのグラフラプラシアンに、ハイパーグラフラプラシアンは1つのハイパーグラフに、MC-MI-HOOIは複数のグラフラプラシアン(3階テンソル)に対して適用される。

これまでは次数という言葉をあるノードに接続されたエッジの本数(または重みの総和)として使っていた。一方、あるエッジが幾つのノードと繋がっているかはエッジの次数として区別される。エッジejの次数d(ej)は、接続行列Hgj列目の和として求められる。

  

(41) d ( e j ) = i | V | H g ( i , j ) = 2

ただし、意味的にどのエッジも次数が2なのは自明である。同様に、ハイパーグラフのハイパーエッジhjの次数dhj)は接続行列Hhyperj列目の和として求められる。

  

(42) d ( h j ) = i | V | H hyper ( i , j )

ただし、通常のグラフと異なり、ハイパーエッジの次数は2とは限らない。

式(33)と式(36)から、グラフラプラシアンは、接続行列Hgを用いて、以下のように書き換えることができる。

  

(43) L = 2 D v , weight H g W H g T

ここで、さらにHgHhyperに置き換えたものを、ハイパーグラフラプラシアンという[75, 76, 77, 78]。

  

(44) = 2 D v , weight H hyper W H hyper T

通常Wの対角要素を前述のエッジ次数の逆数 D e 1 をかけることでハイパーエッジごとの正規化を行い、さらに(対称)正規化でノードごとの正規化を行った、以下の式が利用されることが多い。

  

(45) sym = I 1 2 D v 1 2 H hyper D e 1 H hyper T D v 1 2

このハイパーグラフラプラシアンに固有値分解を適用し、固有ベクトルを次元圧縮に利用した場合は、ハイパーグラフラプラシアン固有マップ、さらにそれら固有ベクトルに対してクラスタリング手法を適用したものは、ハイパーグラフスペクトラルクラスタリングと言う[75, 76, 77, 78]。

バイオインフォマティクス分野でハイパーグラフスペクトラルクラスタリングが利用された例としては、Hi-Cの類似実験技術であるPore-Cのデータ解析[79]や、通常の行列分解に正則化項としてハイパーグラフラプラシアンが加わったもので、遺伝子発現データの解析に活用されたものがある[80]。

ハイパーグラフラプラシアンは、多対多の関係性や、状況に応じてパートナーが切り替わるようなグラフを表現しやすいが、CAMCAの関係性と同様、最終的には行列状のLLsymに情報を押しつぶしてしまっており、個々のデータの個体差を考慮したモデリングとはいえない。個体ごとに異なるグラフラプラシアンがあった場合に、それらをまとめて1つのテンソルにして、テンソル分解を適用した例があり[81]、INDSCALの場合と同様、1つのハイパーグラフラプラシアンに平均化した場合と比較して、解析結果が向上する可能性がある。

おわりに

今回は、これまで紹介した計測データの行列・テンソル分解からさらに発展して、分割表や、距離行列、グラフラプラシアンに対しての行列・テンソル分解を紹介した。分割表における行列分解であるCAや、距離行列における行列分解であるMDS、グラフラプラシアンにおける行列分解であるラプラシアン固有マップ、スペクトラルクラスタリングは、既にバイオインフォマティクスの様々な問題に適用されている。今後、これら手法を行列同時分解や、テンソル分解に拡張することによって、同時に考慮する変数を増やし、変数ごとの属性の違いや、個体差を考慮したモデルを構築できるため、それにより解析結果が向上する可能性がある。

これまでの連載で紹介した行列・テンソル分解手法は、確率モデルとして解釈した場合、全てのデータが互いに独立に同じ分布から生成されるindependent and identically distributed(i.i.d.)という仮定に基づいている。i.i.d.は目的関数を簡略化でき最適化しやすいメリットがあるため、多くのモデルが暗にこれを仮定している。一方で、時系列データや空間データでは、時間や空間座標の変化に対して、観測値が連続的に変化するため、隣り合うデータは似た値になりやすい。このようなデータにおいては、i.i.d.の仮定は不適切で、モデルの性能が落ちる可能性があることから、隣接するデータからの影響もモデルに組み込む必要が出てくる。次回は、このような系列データに特化した行列・テンソル分解手法について紹介する。

References
略語リスト

・SNP

Single Nucleotide Polymorphism

・PCA

Principal Component Analysis

・CCA

Canonical Correlation Analysis

・MDS

Multi-Dimensional Scaling

・CA

Correspondence Analysis

・SVD

Singular Value Decomposition

・CMH検定

Cochran-Mantel-Haenszel検定

・MCA

Multiple Correspondence Analysis

・KL距離

Kullback-Leibler Divergence

・PCoA

Principal Coordinate Analysis

・ICA

Independent Component Analysis

・NMF

Non-negative Matrix Factorization

・INDSCAL

Individual Differences Scaling

・ALS

Alternating Least Squares

・i.i.d.

independent and identically distributed

著者略歴

露崎 弘毅
2015年、東京理科大学生命創薬科学科博士後期課程終了。博士(薬科学)。同年より、理化学研究所情報基盤センターバイオインフォマティクス研究開発ユニット(現在の所属1)に在籍。現在、基礎科学特別研究員及びJSTさきがけ研究員を兼任。パッケージングに特化したハッカソンBio “Pack” athonを主催。
ホームページ:https://sites.google.com/view/kokitsuyuzaki

 
© 2022 日本バイオインフォマティクス学会

This article is licensed under a Creative Commons [Attribution-NonCommercial-ShareAlike 4.0 International] license.
https://creativecommons.org/licenses/by-nc-sa/4.0/
feedback
Top