Rによる階層的クラスタリング
階層的クラスタリングを行う．以下のような csv ファイルにあるデータをクラスタリングする．以下は20個のデータからなり，それぞれが f1 から f10 の10個の値を持っている．クラスタ解析をすることで，これらの20個のデータのどれとどれが近いか，遠いかがわかる．データは20種類のアミノ酸に対する値．
,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10 ALA,-1.56,-1.67,-0.97,-0.27,-0.93,-0.78,-0.2,-0.08,0.21,-0.48 GLU,-1.45,0.19,-1.61,1.17,-1.31,0.4,0.04,0.38,-0.35,-0.12 MET,-1.4,0.18,-0.42,-0.73,2,1.52,0.26,0.11,-1.27,0.27 LEU,-1.04,0,-0.24,-1.1,-0.55,-2.05,0.96,-0.76,0.45,0.93 VAL,-0.74,-0.71,2.04,-0.4,0.5,-0.81,-1.07,0.06,-0.46,0.65 ILE,-0.73,-0.16,1.79,-0.77,-0.54,0.03,-0.83,0.51,0.66,-1.78 GLN,-0.47,0.24,0.07,1.1,1.1,0.59,0.84,-0.71,-0.03,-2.33 HIS,-0.41,0.52,-0.28,0.28,1.61,1.01,-1.85,0.47,1.13,1.63 LYS,-0.34,0.82,-0.23,1.7,1.54,-1.62,1.15,-0.08,-0.48,0.6 PHE,-0.21,0.98,-0.36,-1.43,0.22,-0.81,0.67,1.1,1.71,-0.44 CYS,0.12,-0.89,0.45,-1.05,-0.71,2.41,1.52,-0.69,1.13,1.1 ARG,0.22,1.27,1.37,1.87,-1.7,0.46,0.92,-0.39,0.23,0.93 THR,0.26,-0.7,1.21,0.63,-0.1,0.21,0.24,-1.15,-0.56,0.19 TRP,0.3,2.1,-0.72,-1.57,-1.16,0.57,-0.48,-0.4,-2.3,-0.6 ASP,0.58,-0.22,-1.58,0.81,-0.92,0.15,-1.52,0.47,0.76,0.7 SER,0.81,-1.08,0.16,0.42,-0.21,-0.43,-1.89,-1.15,-0.97,-0.23 ASN,1.14,-0.07,-0.12,0.81,0.18,0.37,-0.09,1.23,1.1,-1.73 TYR,1.38,1.48,0.8,-0.56,0,-0.68,-0.31,1.03,-0.05,0.53 GLY,1.46,-1.96,-0.23,-0.16,0.1,-0.11,1.32,2.36,-1.66,0.46 PRO,2.06,-0.33,-1.15,-0.75,0.88,-0.45,0.3,-2.3,0.74,-0.28
最初に，このデータをコピーし，以下のコマンドで変数 dx に読み込む．エクセルからもコピーして読み込める．その場合は，sep="," は指定しない．
$ dx=data.frame(read.table(file="clipboard",sep=",",header=T,row.names=1))
または以下のように，ファイルから読み込む．ここで，x.csv はファイル名．コマンドはファイルが存在しているディレクトリで打つ必要がある．
$ dx=data.frame(read.table(file="x.csv",sep=",",header=T,row.names=1))
このように生成した変数 dx に対して，以下のコマンドを打つことで各変数間の距離を計算するために，以下のコマンドを打つ．この rd には，各データ間の距離が格納される．
$ rd=dist(dx)
この距離を用いて階層的クラスタリングを行うに以下のコマンドを打つ．ここで，method にはクラスタリングする際のアルゴリズムの名前を指定する．最も一般的なのは average とか ward.D2 となる．
$ rc=hclust(d=rd,method="ward.D2")
コマンド hclust() で利用することができるアルゴリズムには以下のようなものがある．生物学分野では群平均法が最も用いられてきた．完全連結法はより左右対称のバランス様の樹形図を生成するし，単連結法はより鎖状の樹形図を生成する．ウォード法は最も外れ値の影響を受けにくい．様々なアルゴリズムが使用できるが，クラスタリングはデータの縮約を目的にした発見的または探索的な手法であり，データによって使われるべきアルゴリズムが決まっているわけではない．すなわち，クラスタリングの結果からは一定の示唆は与えるモノであるが，それ自体が何かを確定的に主張するものではない．なので，実際の計算の際は色々な評価基準を試してみるべきである．
| 評価基準 | 詳細 |
|---|---|
| single | 単連結法 (最短距離法)：クラスター間で各クラスターに属するOTU間の距離が最短のものをそのクラスターの距離とする． |
| complete | 完全連結法 (最長距離法)：クラスター間で各クラスターに属するOTU間の距離が最長のものをそのクラスターの距離とする． |
| average | 群平均法：UPGMA とも言う．各クラスターに属する全OTU間の距離の平均をクラスター間の距離とする． |
| centroid | 重心法：クラスター間の重心間距離からクラスター間距離を導く．別のクラスターとの距離は重心間距離を要素数で重み付けした点から生成する． |
| median | メディアン法：重心法と若干異なり，重心間距離を要素数でウエイトせずに単純に中点から別のクラスターとの距離を生成する． |
| mcquitty | McQuitty 法：クラスターAとクラスターBを併合した新クラスターCと別のクラスターDの距離を距離ADおよびBDから決定する． |
| ward.D | Ward 法 (最小分散法)：クラスター内の OTU 間の距離平方和の増加量が最小になるクラスターを併合させる．R のバージョン3.0.3以前は単に ward のみで指定した． |
| ward.D2 | Ward 法 (最小分散法)：クラスター内の OTU 間の距離平方和の増加量が最小になるクラスターを併合させる．ward.D は本来の Ward 法の実装ではなく，こっちが本物．ward.D では本来すべき距離の二乗が計算されない． |
上で得られたクラスタリングの結果からデンドログラムを生成するには以下のように打つ．
$ plot(rc)
これを実行した結果は以下のようになる．
PHE と TYR，THR と SER，VAL と ILE が近くなっているのは実際のアミノ酸の性質を考えても納得できる．このように生成したクラスタリング結果に対して，以下のコマンドを打つと，指定したクラスタ数で分割した場合のグループのエントリー (値) を知ることができる．
$ cutree(tree=rc,k=2)
これを実行した結果は以下のようになる．ALA，GLU，MET，LEU，HIS，LYS，PHE，TRP，ASP，TYR，GLY の11個がクラスタ1，残りがクラスタ2に分類される．
$ ALA GLU MET LEU VAL ILE GLN HIS LYS PHE CYS ARG THR TRP ASP SER ASN TYR GLY PRO
1 1 1 1 2 2 2 1 1 1 2 2 2 1 1 2 2 1 1 2
また，以下のコマンドを打つと，デンドログラムの高さを分類の閾値に設定することができる．
$ cutree(tree=rc,h=5)
これを実行した結果は以下のようになる．
$ ALA GLU MET LEU VAL ILE GLN HIS LYS PHE CYS ARG THR TRP ASP SER ASN TYR GLY PRO
1 1 2 3 4 4 5 1 3 3 6 6 4 2 1 4 5 3 7 5
リファレンス
- Rackovsky S and Scheraga HA, On the Information Content of Protein Sequences, Journal of Biomolecular Structure & Dynamics, 28(4):593-594, 2011