ホームへ
script03.nb.html
File Size: 927 kb
File Type: html
ファイルのダウンロード
教師なし学習（クラスタ－分析） 
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))
}

LS0tDQp0aXRsZTogIouzjnSCyIK1inePS4Fpg06DiYNYg16BfJWqkM2BaiINCmF1dGhvcjogIll1c3VrZSBNYXRzdWkiDQphZmZpbGlhdGlvbjogIpa8jMOJrpHlineDVoNYg2WDgJC2laiKd5WqluwiDQpkYXRlOiAijcWPSY1YkFaT+oFGMjAxN5RONYyOMTeT+iINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCmBgYHtyLGV2YWw9Rn0NCmhvbWUgPC0gW5R6lXqCtYK9g3SDSIOLg1+C8INSg3OBW4K1gr2DcINYXQ0Kc2V0d2QoaG9tZSkNCmBgYA0KYGBge3IsZWNobz1GfQ0KbGlicmFyeShrbml0cikNCmBgYA0KDQqDToOJg1iDXoFblaqQzYLGgs2BQYFpjpaRT4LMj+6V8YLwjmeC7YK4gsmBapWhkJSCzIzCkcyC8JWql96BaYNPg4uBW4Nzg5ODT4FqgreC6ZX7lkCCxYK3gUKI4pNgjnGUrYy7g2aBW4NegsyX4YLwjWyCpoLpgsaBQZStjLuDcINegVuDk4LMjpeCxIKigumI4pNgjnGMUYLwk6+S6IK1gr2CooFBgqCC6YKigs2DVIOTg3aDi4xRgvCTr5LogrWCvYKij+qNh4LIgseCyZNLl3CCt4LpgrGCxoKqgsWCq4LcgreBQoNOg4mDWINegVuVqpDNgsWCzYzCkcyK1ILMYJfejpeTeGCC8I1sgqaCxIFBl96Ol5N4gsyNgoKijMKRzJOvjm2CzZOvgraDT4OLgVuDdoFpg06DiYNYg16BW4FqgsmCyILpguaCpILJgrWCxIFBg06DiYNYg16BW4LwjGCQrIK1gtyCt4FCl+GCxoK1gsSBQYNUg5ODdoOLkJSCqoJRgsWBQTMwjMKCzIjik2COcYLMlK2Mu4NmgVuDXoLwjWyCpoLEgt2C3IK1guWCpIFCDQoNCmBgYHtyfQ0Kc2V0LnNlZWQoMTAwKQ0KbiA8LSAxMA0KcCA8LSAyDQptdSA8LSAxMDANCnNpemUgPC0gMTAwDQp4IDwtIG1hdHJpeChybmJpbm9tKG4qcCxzaXplID0gc2l6ZSxtdSA9IG11KSxuY29sPXApDQp5IDwtIG1hdHJpeChybmJpbm9tKG4qcCxzaXplID0gc2l6ZSxtdSA9IG11ICsgMjAwKSxuY29sPXApDQp6IDwtIG1hdHJpeChybmJpbm9tKG4qcCxzaXplID0gc2l6ZSxtdSA9IG11ICsgNDAwKSxuY29sPXApDQpkYXRhX2V4YW1wbGUgPC0gcmJpbmQoeCx5LHopDQpjb2xuYW1lcyhkYXRhX2V4YW1wbGUpIDwtIHBhc3RlMCgiU2FtcGxlIiwxOnApDQpyb3duYW1lcyhkYXRhX2V4YW1wbGUpIDwtIHBhc3RlMCgiR2VuZSIsc2VxKDEsbiozKSkNCmthYmxlKGRhdGFfZXhhbXBsZVsxOjEwLF0pDQpwbG90KGRhdGFfZXhhbXBsZSx4bGFiID0gIlNhbXBsZSAxIix5bGFiPSJTYW1wbGUgMiIsbWFpbj0iR2VuZSBleHByZXNzaW9uIG9mIHNhbXBsZSAxIGFuZCAyIikNCmBgYA0KjKmCvZbagsmCzYFBglKCwoLMg06DiYNYg16BW4KqgqCC6IK7gqSCxYK3gUKOwI3bgsmDToOJg1iDXoFblaqQzYLwjsCNc4K1gsSC3YLcgrWC5YKkgUKCoILGgsWBQZDglr6Ct4Lpk/GCwoLMlfuWQILwgsaC6IKggqaCuI7AjXOCtYLEgt2C3IK1guWCpIFCDQoNCmBgYHtyfQ0KI5dcgt+DToOJg1iDXoFbkJSC8I53kuiCt4LpDQprIDwtIDMNCg0KI4NOg4mDWINegVuVqpDNKIpLkXeMXoNOg4mDWINegVuVqpDNlkApDQpkIDwtIGRpc3QoZGF0YV9leGFtcGxlKQ0KcmVzdWx0X2hjbHVzdCA8LSBoY2x1c3QoZCkNCnJlc3VsdF9jbHVzdGVyIDwtIGN1dHJlZShyZXN1bHRfaGNsdXN0LGspDQpgYGANCg0KjIuJyoLwg3aDjYNig2eCtYLEgt2C3IK3gUINCg0KYGBge3J9DQpjb2xfa2V5cyA8LSBjKCJyZWQiLCJibHVlIiwicHVycGxlIikNCg0KY29sX2xhYmVsIDwtIGNvbF9rZXlzW3Jlc3VsdF9jbHVzdGVyXQ0KcGxvdChkYXRhX2V4YW1wbGUsY29sPWNvbF9sYWJlbCx4bGFiID0gIlNhbXBsZSAxIix5bGFiPSJTYW1wbGUgMiIsbWFpbj0iUmVzdWx0IG9mIGNsdXN0ZXJpbmciKQ0KbGVnZW5kKCJib3R0b21yaWdodCIsbGVnZW5kPXBhc3RlKCJDbHVzdGVyIiwxOmspLHBjaD1yZXAoMSxrKSxjb2w9Y29sX2tleXMpDQpgYGANCg0KjIuJyoLNjcWPiYLMl1yReoLMksqC6ILJglKCwoLMg06DiYNYg16BW4LJlaqCr4LpgrGCxoKqgsWCq4LcgrWCvYFCg06DiYNYg16BW5CUgvCCoILngqmCtoLfjIiC34LIgq+C6oLOgsiC54LIgqKVlJWqgs2PrYK1i0OCyYLIguiC3IK3gqqBQYK7guqCzYzjgsWLY5hfgrWC3IK3gUINCg0Kg06DiYNYg16BW5WqkM2WQILNikuRd4xegsaU8YpLkXeMXoLMk/GCwoLJkeWVyoLFgquC3IK3gUKKS5F3jF6CzYFBgqCC6YzCkcyCyY3FguCX3o6XgrWCvZXKgsyMwpHMgvCSVIK1gsSBQY6fgViCxoNPg4uBW4Nzg5ODT4LwgrWCxIKigq2V+5ZAgsWCt4FClPGKS5F3jF6CzYFBl1yC34yIgt+CvZCUgsyR45Vck0mCyJNfgsmRzoK1gsSBQY3FguCL34KijMKRzILwkFWC6JWqgq+CxIKigq2V+5ZAgsWCt4FCDQoNCopLkXeMXoNOg4mDWINegVuVqpDNgsWCzYFBgtyCuIzCkcyK1ILMl96Ol5CrgWmOwI3bgs2QlIp3k0mCyYi1gqKC4oK3gqKC5oKkgsmU8ZfejpeQq4KggumCooLNi5eXo4FqgvCMdo5agrWC3IK3gUKCu4LMgr2C34LJgUFSgsWCzWBkaXN0YIrWkJSCqpdwiNOCs4LqgsSCooLcgreBQg0KYGBge3J9DQpkIDwtIGRpc3QoZGF0YV9leGFtcGxlKQ0KYGBgDQpgYGB7cn0NCmRfZGlzcCA8LSBhcy5tYXRyaXgoZCkNCmthYmxlKGRfZGlzcFsxOjEwLDE6MTBdKQ0KYGBgDQoNCo6fgsmBQYKxgsyLl5ejgvBgaGNsdXN0YIrWkJSCyZP8l82CtYLcgreBQg0KDQpgYGB7cn0NCihyZXN1bHRfaGNsdXN0IDwtIGhjbHVzdChkKSkNCmBgYA0KikuRd4xeg06DiYNYg16BW5WqkM2CzIyLicqC8IypgumR5ojqleCCzYFBjveP85B9gWmDZoOTg2iDjYNPg4mDgDtEZW5kcm9ncmFtgWqCxozEgs6C6oLpkH2C8IrPjkCCt4LpgrGCxoLFgreBQg0KYGBge3J9DQpwbG90KHJlc3VsdF9oY2x1c3QpDQpgYGANCg0KgrGCzJB9gs2BQYLHgsyC5oKkgsiPh5TUgsWMwpHMgsyDeYNBgqqDT4OLgVuDc4OTg0+Cs4Lqgr2CqYLwjqaCtYLEgqiC6IFBjYKCs4LNg0+Di4Fbg3ODk4NPgsyN24LJl3CCooLnguqCvZTxl96Ol5N4gsWCt4FCl+GCpoLOgUGJRZJbgsxHZW5lMTGCxkdlbmUxMoLMg3mDQYKqg0+Di4Fbg3ODk4NPgrOC6oLEgUGOn4LJgruCzINPg4uBW4N2gsaNxYLgi9+CooLMgs1HZW5lMTmCvoLBgr2CsYLGgvCVXIK1gsSCooLcgreBQoNmg5ODaIONg0+DiYOAgs2DToOJg1iDXoFbgsyNXJGigvCVXIK1gsSCooLcgreBQoLCgtyC6IFBgseCsYLFkIWVvYLJg2aDk4Nog42DT4OJg4CC8GCQ2JJmYIK3guqCzoFBYONZl+2CyGCVqoKpguqV+4LwgreC6YKpgvCOpoK1gsSCooLcgreBQo2hifGCzYJSgsKCzINPg4uBW4N2gqqCoILogruCpILIgsyCqoLtgqmC6ILcgreBQoLIgrqCyILngUFHZW5lMSAtIEdlbmUxMILMg0+Di4Fbg3aCxkdlbmUxMSAtIEdlbmUyMILMg0+Di4Fbg3aCqoNPg4uBW4Nzg5ODT4KzguqC6YLGgquCzI2CgrOCzYFBimWDT4OLgVuDdpPggsyI4pNgjnGK1ILMlPGX3o6XkKuC5oLoguCU8Y/tgsmR5YKrgqKCrYLIgsGCxIKoguiBQYLcgr1HZW5lMSAtIEdlbmUyMILMg0+Di4Fbg3aCxkdlbmUyMSAtIEdlbmUzMILMg0+Di4Fbg3aK1ILMlPGX3o6XkKuC4JOvl2yCvoKpgueCxYK3gUINCmBgYHtyLGVjaG89Rn0NCnBsb3QocmVzdWx0X2hjbHVzdCkNCmFibGluZShoPTIwMCxjb2w9InJlZCIpDQpgYGANCg0KDQqCsYLMg2aDk4Nog42DT4OJg4CC8JdwgqKCxIFBglKCwoLMg06DiYNYg16BW4LJlaqCr4LpgsmCzYFBYGN1dHJlZWCK1pCUgvCOZ4KigtyCt4FCDQpgYGB7cn0NCmNsdXN0ZXJzIDwtIGN1dHJlZShyZXN1bHRfaGNsdXN0LGspDQpgYGANCg0KYGBgYHtyLGVjaG89Rn0NCmthYmxlKG1hdHJpeChjbHVzdGVycyxucm93PTEsZGltbmFtZXM9bGlzdCgiQ2x1c3RlciBMYWJlbCIscGFzdGUoIkdlbmUiLHNlcSgxLG4qMykpKSkpDQpgYGANCg0KikuRd4xeg06DiYNYg16BW4nwkM2CyYLNgUGCooKtgsKCqYLMg2+DioNHgVuDVoOHg5OCqoKgguiC3IK3gUINCmBgYHtyLGZpZy53aWR0aD0yMH0NCnJlc3VsdF9oY2x1c3QxIDwtIGhjbHVzdChkLG1ldGhvZD0ic2luZ2xlIikNCnJlc3VsdF9oY2x1c3QyIDwtIGhjbHVzdChkLG1ldGhvZD0iY29tcGxldGUiKQ0KcmVzdWx0X2hjbHVzdDMgPC0gaGNsdXN0KGQsbWV0aG9kPSJhdmVyYWdlIikNCnJlc3VsdF9oY2x1c3Q0IDwtIGhjbHVzdChkLG1ldGhvZD0id2FyZC5EIikNCnBhcihtZnJvdz1jKDIsMikpDQpwbG90KHJlc3VsdF9oY2x1c3QxKQ0KcGxvdChyZXN1bHRfaGNsdXN0MikNCnBsb3QocmVzdWx0X2hjbHVzdDMpDQpwbG90KHJlc3VsdF9oY2x1c3Q0KQ0KYGBgDQoNCg0KjMKRzILMg3mDQYLwg0+Di4Fbg3ODk4NPgrWCxIFBjp+CyYK7gsyDT4OLgVuDdoLGjcWC4IvfgqKMwpHMgWmCoILpgqKCzYNPg4uBW4N2gWqC8JJUgreCxoKrgsmCzYFBiOqCwoLMg0+Di4Fbg3aCyYLNlaGQlILMjMKRzIKqityC3ILqgsSCooLpgsyCxYFBg0+Di4Fbg3aCxoNPg4uBW4N2gsyU8ZfejpeTeILMjHaOWoLMjmSV+4LwjIiC34LIgq+C6oLOgsiC6ILcgrmC8YFCl+GCpoLOgUGI6oLCltqCzGBtZXRob2Q9InNpbmdsZSJggs2NxZJai5eXo5ZAgsaMxILOguqBQYNPg4uBW4N2itSCyYKogq+C6YK3gteCxILMjMKRzIrUgsyU8ZfejpeTeILMgqSCv43Fj6yCzILggsyC8INPg4uBW4N2itSCzIuXl6OCxoK1gsSXcIKigumV+5ZAgsWCt4FCYG1ldGhvZD0iY29tcGxldGUiYILGYG1ldGhvZD0iYXZlcmFnZSJggs2NxZHli5eXo5ZAgqiC5oLRlb2Lz4uXl6OWQILGjMSCzoLqgUGDT4OLgVuDdorUgsmCqIKvgumCt4LXgsSCzIzCkcyK1ILMlPGX3o6Xk3iCzI3FkeWSbILGlb2Lz5JsgvCXcIKigtyCt4FCYG1ldGhvZD0id2FyZC5EImCCzYNFg0iBW4NolkCCxozEgs6C6oFBg0+Di4Fbg3aK1ILFgsyVqo5VgsyRnZWqgvCU8ZfejpeTeILGgrWCxJdwgqKCxIKigtyCt4FCgUKDRYNIgVuDaJZAgs2BQY/Qie6CtYK9lfuWQILMkoaCxZdCiOqBQYNmgVuDXoLMk52MdpNJgsiCzoLngsKCq4LwjWyXtoK1gr2V+5ZAgsWCt4FCgseCzJX7lkCCqpfHgqKCzIKpgsmCwoKigsSBQYLmgq2Ov5bigvCO84KvgtyCt4KqgUGTmoKmgs1ggu2CqYLngsiComCCxYK3gUKNoYnxgsyX4YLFgs2BQYLHgsyV+5ZAgsWC4IyLicqCzZXPgu2C6ILcgrmC8YFCgsKC3ILogUGWe5OWgsyDToOJg1iDXoFbjVyRooLGgqKCpILMgs2BQZX7lkCCyYjLkbaCtYLEg1KDjYNSg42CxpXPgu2C6YLggsyCxYLNgsiCooKxgsaC8I6mjbSCtYLEgqKC3IK3gUINCmBgYHtyLGVjaG89Rn0NCmthYmxlKHJiaW5kKHNpbmdsZT1jdXRyZWUocmVzdWx0X2hjbHVzdDEsayksY29tcGxldGU9Y3V0cmVlKHJlc3VsdF9oY2x1c3QyLGspLGF2ZXJhZ2U9Y3V0cmVlKHJlc3VsdF9oY2x1c3QzLGspLHdhcmQuRD1jdXRyZWUocmVzdWx0X2hjbHVzdDQsaykpKQ0KYGBgDQoNCoK1gqmCtYFBguCCv4LrgvGCu4LqgryC6oKqlrOI05ahgsmS8YjEgrOC6oK9gu2Cr4LFgs2CyIKtk8GSpYLNjp2CwYLEgqKC3IK3gUKNxZJai5eXo5ZAgs2SUI+DgsmL34KiguCCzILwjp+BWILGgtyCxoLfgsSCooKtgsyCxYFBkVORzILGgrWCxILNgUGL34KiguCCzI+HgsmVwILUjFiM/IKqgqCC6ILcgreBQpPBgsmWvoptgsiDT4OLgVuDdo1ckaKCqoLIgqKCxoKrgsmBQYjqgsKCzINPg4uBW4N2gsmI6oLCgsyMwpHMgqqOn4FYgsmDT4OLgVuDc4OTg0+Cs4LqgsSI6oLCgsyCqYK9gtyC6ILJgsiCwYLEgrWC3IKkjLuP24FplaqX3pedmF+CxYLNg2CDRoFbg5OMu4/bgsaMxILOguqC6YFqgqqLToKxguiC3IK3gUKI6pX7gUGNxZK3i5eXo5ZAgs2DYINGgVuDk4y7j9uC8InxlPCCtYLEgUGI6oLCgsyDT4OLgVuDdoLJgqmCvYLcgueCyIKiguaCpILJg06DiYNYg16BW4LwjGCQrIK3gumMWIz8gqqCoILogtyCt4FCg0WDSIFbg2iWQILNguCCxoLggsaCzINmgVuDXoLOgueCwoKrgsyNXJGigvCUvYlmgrWCvYxggsWDToOJg1iDXoFbjGCQrIK3gumCzILFgUGI6pTUjcWPiYLJjKmCvZB9gsmSvIrPk0mCyYvfgqKMYILFg0+Di4Fbg3ODk4NPgreC6YxYjPyCqoKgguiC3IK3gUINCg0KDQoNCoKzgsSBQYNOg4mDWINegVuQlILJgsKCooLEgs2BQYKxguqC3ILFl1yC35degqaC6YLGgqKCpJFPkvGCxZhigvCQaYLfgsSCq4LcgrWCvYFCjsCN24FBg06DiYNYg16BW5WqkM2CzWCLs450gsiCtYp3j0tggsaCooKkgq2C54KigsWCt4KpgueBQZCzifCCzINOg4mDWINegVuQlILIgseCzYLIgqKCxoKigqSCzIKqkLOCtYKigsaMvoKmgtyCt4FCgqKCrYLCgqmCzINBg3aDjYFbg2CCqpLxiMSCs4LqgsSCooLcgreBQojqgsKCzYKxguqC3ILFi2OYX4K1gr2DToOJg1iDXoFblaqQzZZAgsaCzZFTgq2I2YLIgumDQYN2g42BW4NggsWBQYOCg2aDi4N4gVuDWINOg4mDWINeg4qDk4NPgsaMxILOguqBQYplg06DiYNYg16BW4LGgs2Ud4zjgsmCoILpiNmCyILpleqPV5JjgvCVXIK3gsaNbIKmgsSBQYNmgVuDXoLwjcWC4ILmgq2Q4Ja+greC6YLmgqSCyIptl6aVqpV6gvCTloLEgs2C34LplfuWQILFgreBQoOCg2aDi4N4gVuDWINOg4mDWINeg4qDk4NPgs2DToOJg1iDXoFbkJSC8I6pk66CxYyIgt+C6YKxgsaCqonClFyCxYK3gUKJupB9gs2Ud4zjgsmCooKtgsKCqYLMkLOLS5WqlXqCqoKggumCxom8kuiCtYLEk5aCxILNgt+CvY/qjYeCxYFBjk+CwoLMg06DiYNYg16BW4Kqk6+S6IKzguqCxIKigtyCt4FCDQoNCmBgYHtyfQ0KbGlicmFyeShtY2x1c3QpDQpyZXN1bHRfbWNsdXN0IDwtIE1jbHVzdChkYXRhX2V4YW1wbGUpDQpjb2xfbGFiZWwgPC0gY29sX2tleXNbcmVzdWx0X21jbHVzdCRjbGFzc2lmaWNhdGlvbl0NCnBsb3QocmVzdWx0X21jbHVzdCx3aGF0ID0gImRlbnNpdHkiKQ0KcG9pbnRzKGRhdGFfZXhhbXBsZSxjb2w9Y29sX2xhYmVsKQ0KYGBgDQoNCoLggqSI6oLCgs2BQZO+gueC6oK9g06DiYNYg16BW4LMYIptgqmCs2CC8JLol8qJu4K1gsSVXYm/greC6YNBg3aDjYFbg2CCxYK3gUKNoYnxgs2BQZPxgsKCzJX7lkCC8I/Qie6CtYLcgreBQojqgsKW2oLNYHNpbGhvdWV0dGVggsaMxILOguqC6YLggsyCxYFBjGCQrIKzguqCvYNOg4mDWINegVuCqoLHgsyS9pN4gUGCtYLBgqmC6JWql6OCtYLEgqKC6YLMgqmC8JVdib+Ct4LpjneVV4LFgreBQopliOKTYI5xgrKCxoLJjLuN3YLMg06DiYNYg16BW4LJka6CtYLEgqKC6YptgqmCs4LwLTEojLuN3YLMg06DiYNYg16BW4LFgs2CyIKiKYKpgucxKIy7jd2CzINOg4mDWINegVuCxYKggukpgsySbILFlV2Jv4K1gtyCt4FCiMiJuoLJg06DiYNYg16BW5CUgqoygqmC5zWC3ILFgsxzaWxob2V0dGUgcGxvdILwjqaCtYLcgreBQoplll+CqopliOKTYI5xgsmRzoK3gulzaWxob3VldHRlkmyC8JVcgrWCxIKoguiBQTGCyYvfgq+C6oLOjLuN3YLMg06DiYNYg16BW4LJka6CtYLEgqKC6YnClFyQq4KqjYKCooLGifCO34K1gtyCt4FCimWQfYLMiUWRpILJgqCC6YLMgqqBQYNOg4mDWINegVuCsoLGgsxzaWxob2V0dGWSbILMlb2Lz4LFgUGCsYLqgqqNgoKvguqCzoNOg4mDWINegVuCzIptgqmCs4LNjYKCooLGjWyCpoLcgreBQoLcgr2BQZFTkcyCxoK1gsSCzYFBibqCyYKggulgQXZlcmFnZSBzaWxob3VldHRlIHdpZHRoYIKqkVORzILMc2lob3VldHRlkmyCzJW9i8+CxYFBgrGCzJJsgqqNgoKig06DiYNYg16BW5CUgqqNxY9Jk0mCyWCXx4KiYINOg4mDWINegVuQlIK+gsaNbIKmgtyCt4FCjaGJ8YLNaz0zgsaCyILogtyCt4FCDQoNCoLggr+C64LxgqKCrYLCgqmCzINOg4mDWINegVuQlILGgsyK1ILFlPeWrYLIj+qNh4LJgs2BQYNOg4mDWINegVuCsoLGgsxzaWhvdXR0ZZJsgvCP2o3XgsmMqYLEifCO34K3gumVS5d2guCCoILogtyCt4FCaz00LDWCxYLNlYmCzHNpbGhvdWV0dGWSbILwjp2Cwojik2COcYKqityC3ILqgsSCooK9guiBQXNpaG91ZXR0ZYLMlb2Lz5JsgqqS4YKig06DiYNYg16BW4KqgqCC6YLMgsWBQYKggtyC6JfHgqKCxoy+gqaC3IK5gvGBQoK7grGCxWs9MiwzgsyCx4K/gueCqYLGgsiC6IFBaz0ygsWCzYjqgsKW2oLMg06DiYNYg16BW4LMc2lob3VldHRlkmyCqpLhgq2BQYK7guqC54LwgrOC54LJlaqCr4K9az0zgsWCzY2Cgq2CyILBgsSCooLpgsyCxYFBaz0zgsaJ8I7fguCCxYKrgtyCt4FCDQpgYGB7cixmaWcud2lkdGg9MTAsZmlnLmhlaWdodD0xMH0NCmxpYnJhcnkoY2x1c3RlcikNCmtzIDwtIDI6NQ0KcGFyKG1mcm93PWMoMiwyKSkNCmZvcihpIGluIHNlcV9hbG9uZyhrcykpew0KICBrIDwtIGtzW2ldDQogIGNscyA8LSBjdXRyZWUocmVzdWx0X2hjbHVzdCxrc1tpXSkNCiAgcGxvdChzaWxob3VldHRlKGNscyxkKSxtYWluPXBhc3RlMCgiaz0iLGspKQ0KfQ0KYGBgDQoNCg==
Proudly powered by Weebly