教師なし学習（クラスタ－分析）

home <- [配布したフォルダをコピーしたパス]
setwd(home)

クラスター分析とは、（事前の情報を使わずに）複数の個体を分類（グルーピング）する方法です。遺伝子発現データの例を考えると、発現パターンの似ている遺伝子群を同定したい、あるいはサンプル群を同定したい場合などに適用することができます。クラスター分析では個体間の類似度を考えて、類似度の高い個体同士は同じグループ（クラスター）になるようにして、クラスターを形成します。例として、サンプル数が２で、30個の遺伝子の発現データを考えてみましょう。

set.seed(100)
n <- 10
p <- 2
mu <- 100
size <- 100
x <- matrix(rnbinom(n*p,size = size,mu = mu),ncol=p)
y <- matrix(rnbinom(n*p,size = size,mu = mu + 200),ncol=p)
z <- matrix(rnbinom(n*p,size = size,mu = mu + 400),ncol=p)
data_example <- rbind(x,y,z)
colnames(data_example) <- paste0("Sample",1:p)
rownames(data_example) <- paste0("Gene",seq(1,n*3))
kable(data_example[1:10,])

	Sample1	Sample2
Gene1	79	83
Gene2	109	90
Gene3	96	108
Gene4	100	96
Gene5	101	100
Gene6	109	114
Gene7	89	108
Gene8	115	95
Gene9	95	105
Gene10	109	93

plot(data_example,xlab = "Sample 1",ylab="Sample 2",main="Gene expression of sample 1 and 2")

見た目には、３つのクラスターがありそうです。実際にクラスター分析を実行してみましょう。あとで、説明する二つの方法をとりあえず実行してみましょう。

In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

#クラスター分析(階層型クラスター分析法)
d <- dist(data_example)
result_hclust <- hclust(d)
result_cluster <- cutree(result_hclust,k)

結果をプロットしてみます。

col_keys <- c("red","blue","purple")
col_label <- col_keys[result_cluster]
plot(data_example,col=col_label,xlab = "Sample 1",ylab="Sample 2",main="Result of clustering")
legend("bottomright",legend=paste("Cluster",1:k),pch=rep(1,k),col=col_keys)

結果は最初の予想の通りに３つのクラスターに分けることができました。クラスター数をあらかじめ決めなければならない部分は少し気になりますが、それは後で議論します。

クラスター分析法は階層型と非階層型の二つに大別できます。階層型は、ある個体に最も類似した別の個体を探して、次々とグルーピングをしていく方法です。非階層型は、予め決めた数の代表的な点に対して、最も近い個体を振り分けていく方法です。

階層型クラスター分析では、まず個体間の類似性（実際は数学的に扱いやすいように非類似性あるいは距離）を計算します。そのために、Rではdist関数が用意されています。

d <- dist(data_example)

d_disp <- as.matrix(d)
kable(d_disp[1:10,1:10])

	Gene1	Gene2	Gene3	Gene4	Gene5	Gene6	Gene7	Gene8	Gene9	Gene10
Gene1	0.00000	30.80584	30.232433	24.698178	27.802878	43.13931	26.925824	37.947332	27.202941	31.622777
Gene2	30.80584	0.00000	22.203603	10.816654	12.806249	24.00000	26.907248	7.810250	20.518285	3.000000
Gene3	30.23243	22.20360	0.000000	12.649111	9.433981	14.31782	7.000000	23.021729	3.162278	19.849433
Gene4	24.69818	10.81665	12.649111	0.000000	4.123106	20.12461	16.278821	15.033296	10.295630	9.486833
Gene5	27.80288	12.80625	9.433981	4.123106	0.000000	16.12452	14.422205	14.866069	7.810250	10.630146
Gene6	43.13931	24.00000	14.317821	20.124612	16.124516	0.00000	20.880613	19.924859	16.643317	21.000000
Gene7	26.92582	26.90725	7.000000	16.278821	14.422205	20.88061	0.000000	29.068884	6.708204	25.000000
Gene8	37.94733	7.81025	23.021729	15.033296	14.866069	19.92486	29.068884	0.000000	22.360680	6.324555
Gene9	27.20294	20.51828	3.162278	10.295630	7.810250	16.64332	6.708204	22.360680	0.000000	18.439089
Gene10	31.62278	3.00000	19.849433	9.486833	10.630146	21.00000	25.000000	6.324555	18.439089	0.000000

次に、この距離をhclust関数に入力します。

(result_hclust <- hclust(d))

Call:
hclust(d = d)

Cluster method   : complete 
Distance         : euclidean 
Number of objects: 30 

階層型クラスター分析の結果を見る第一歩は、樹状図（デンドログラム;Dendrogram）と呼ばれる図を観察することです。

plot(result_hclust)

この図は、どのような順番で個体のペアがグルーピングされたかを示しており、高さはグルーピングの際に用いられた非類似度です。例えば、右端のGene11とGene12のペアがグルーピングされて、次にそのグループと最も近いのはGene19だったことを表しています。デンドログラムはクラスターの構造を表しています。つまり、どこで水平にデンドログラムを切断すれば、綺麗な分かれ方をするかを示しています。今回は３つのグループがありそうなのがわかります。なぜなら、Gene1 - Gene10のグループとGene11 - Gene20のグループがグルーピングされるときの高さは、各グループ内の遺伝子間の非類似性よりも非常に大きいくなっており、またGene1 - Gene20のグループとGene21 - Gene30のグループ間の非類似性も同様だからです。

このデンドログラムを用いて、３つのクラスターに分けるには、cutree関数を使います。

clusters <- cutree(result_hclust,k)

	Gene 1	Gene 2	Gene 3	Gene 4	Gene 5	Gene 6	Gene 7	Gene 8	Gene 9	Gene 10	Gene 11	Gene 12	Gene 13	Gene 14	Gene 15	Gene 16	Gene 17	Gene 18	Gene 19	Gene 20	Gene 21	Gene 22	Gene 23	Gene 24	Gene 25	Gene 26	Gene 27	Gene 28	Gene 29	Gene 30
Cluster Label	1	1	1	1	1	1	1	1	1	1	2	2	2	2	2	2	2	2	2	2	3	3	3	3	3	3	3	3	3	3

階層型クラスター解析には、いくつかのバリエーションがあります。

result_hclust1 <- hclust(d,method="single")
result_hclust2 <- hclust(d,method="complete")
result_hclust3 <- hclust(d,method="average")
result_hclust4 <- hclust(d,method="ward.D")
par(mfrow=c(2,2))
plot(result_hclust1)
plot(result_hclust2)

plot(result_hclust3)
plot(result_hclust4)

個体のペアをグルーピングして、次にそのグループと最も近い個体（あるいはグループ）を探すときには、一つのグループには複数の個体が含まれているので、グループとグループの非類似度の計算の仕方を決めなければなりません。例えば、一つ目のmethod="single"は最短距離法と呼ばれ、グループ間におけるすべての個体間の非類似度のうち最小のものをグループ間の距離として用いる方法です。method="complete"とmethod="average"は最大距離法および平均距離法と呼ばれ、グループ間におけるすべての個体間の非類似度の最大値と平均値を用います。method="ward.D"はウォード法と呼ばれ、グループ間での分散の増分を非類似度として用いています。。ウォード法は、紹介した方法の中で唯一、データの統計的なばらつきを考慮した方法です。どの方法が良いのかについて、よく質問を受けますが、答えはわからないです。今回の例では、どの方法でも結果は変わりません。つまり、本当のクラスター構造というのは、方法に依存してコロコロと変わるものではないことを示唆しています。

	Gene1	Gene2	Gene3	Gene4	Gene5	Gene6	Gene7	Gene8	Gene9	Gene10	Gene11	Gene12	Gene13	Gene14	Gene15	Gene16	Gene17	Gene18	Gene19	Gene20	Gene21	Gene22	Gene23	Gene24	Gene25	Gene26	Gene27	Gene28	Gene29	Gene30
single	1	1	1	1	1	1	1	1	1	1	2	2	2	2	2	2	2	2	2	2	3	3	3	3	3	3	3	3	3	3
complete	1	1	1	1	1	1	1	1	1	1	2	2	2	2	2	2	2	2	2	2	3	3	3	3	3	3	3	3	3	3
average	1	1	1	1	1	1	1	1	1	1	2	2	2	2	2	2	2	2	2	2	3	3	3	3	3	3	3	3	3	3
ward.D	1	1	1	1	1	1	1	1	1	1	2	2	2	2	2	2	2	2	2	2	3	3	3	3	3	3	3	3	3	3

しかし、もちろんそれぞれが無意味に提案されたわけではなく特徴は持っています。最短距離法は単純に近いものを次々とまとめていくので、全体としては、近いもの順に並ぶ傾向があります。特に明確なグループ構造がないときに、一つのグループに一つの個体が次々にグルーピングされて一つのかたまりになってしまう現象（分類理論ではチェーン現象と呼ばれる）が起こります。一方、最長距離法はチェーン現象を回避して、一つのグループにかたまらないようにクラスターを形成する傾向があります。ウォード法はもともとのデータばらつきの構造を反映した形でクラスター形成するので、一番最初に見た図に直観的に近い形でグルーピングする傾向があります。

さて、クラスター数については、これまで予め与えるという前提で話を進めてきました。実際、クラスター分析は教師なし学習というくらいですから、正解のクラスター数などはないというのが正しいと言えます。いくつかのアプローチが提案されています。一つはこれまで議論したクラスター分析法とは全く異なるアプローチで、モデルベースクラスタリングと呼ばれ、各クラスターとは背後にある異なる母集団を表すと考えて、データを最もよく説明するような確率分布を当てはめる方法です。モデルベースクラスタリングはクラスター数を自動で決めることが可能です。下図は背後にいくつかの正規分布があると仮定して当てはめた場合で、三つのクラスターが同定されています。

library(mclust)

    __  ___________    __  _____________
   /  |/  / ____/ /   / / / / ___/_  __/
  / /|_/ / /   / /   / / / /\__ \ / /   
 / /  / / /___/ /___/ /_/ /___/ // /    
/_/  /_/\____/_____/\____//____//_/    version 5.2.2
Type 'citation("mclust")' for citing this R package in publications.

result_mclust <- Mclust(data_example)
col_label <- col_keys[result_mclust$classification]
plot(result_mclust,what = "density")
points(data_example,col=col_label)

もう一つは、得られたクラスターの確かさを定量化して評価するアプローチです。今回は、二つの方法を紹介します。一つ目はsilhouetteと呼ばれるもので、形成されたクラスターがどの程度、しっかり分離しているのかを評価する指標です。各遺伝子ごとに現在のクラスターに属している確かさを-1(現在のクラスターではない)から1(現在のクラスターである)の値で評価します。以下にクラスター数が2から5までのsilhoette plotを示します。各棒が各遺伝子に対するsilhouette値を表しており、1に近ければ現在のクラスターに属している可能性が高いと解釈します。各図の右側にあるのが、クラスターごとのsilhoette値の平均で、これが高ければクラスターの確かさは高いと考えます。また、全体としては、下にあるAverage silhouette widthが全体のsihouette値の平均で、この値が高いクラスター数が最終的に良いクラスター数だと考えます。今回はk=3となります。

もちろんいくつかのクラスター数との間で微妙な場合には、クラスターごとのsihoutte値を詳細に見て解釈する必要もあります。k=4,5では負のsilhouette値を持つ遺伝子が含まれていたり、sihouetteの平均値が低いクラスターがあるので、あまり良いと言えません。そこでk=2,3のどちらかとなり、k=2では一つ目のクラスターのsihouette値が低く、それらをさらに分けたk=3では高くなっているので、k=3と解釈もできます。

library(cluster)
ks <- 2:5
par(mfrow=c(2,2))
for(i in seq_along(ks)){
  k <- ks[i]
  cls <- cutree(result_hclust,ks[i])
  plot(silhouette(cls,d),main=paste0("k=",k))
}