K-means クラスタリング

膨大な数の遺伝子の発現パターンを似たパターン同士にクラスタリングする手法．K-means クラスタリングの原理はこちらの記事がわかりやすい．

>TS<-read.csv("all_genes.csv", header=T) #オブジェクトに発現量データのcsvファイルを指定

>head(TS) #データの確認

>names(TS) #列名の確認

>install.packages("dplyr") #mutate関数（変数変換）のためのパッケージインストール

>library(dplyr)

>TS_log2<-mutate(TS, column1= log2(column1), column2= log2(column2), column3= log2(column3), column4=log2(column4))

#変数のlog2変換（発現ポイントが4つの例）

>head(TS_log2)

>TS_log2[is.nan(TS_log2)]<- -50 #ゼロ発現はlog変換では計算不可能な非数値NaNになるので，とりあえず-50に置換（本当は各値に0.1を足すべきだけどやり方を調査中）

>TS_log2[is.infinite(TS_log2)]<- -50 #-Infも出現，上と同様

>TS_log2<-as.matrix(TS_log2) #class関数でlistと出て変数変換ができない場合は、classをmatrixに変えてやる

>gap<-clusGap(TS_log2[1:1000,],kmeans,K.max = 20,B=100) #最適なクラスター数を機械的に導くためのGap統計量算出．全遺伝子数を扱うのが現実的でない場合は，エクセルのフィルター等で上位発現遺伝子に並び替えて選抜する（この場合だと上位1000遺伝子[1:1000,]，発現量リストを見て恣意的に決定）．全遺伝子でも不可能ではないが計算に一晩くらいかかる

>plot(gap) #最適クラスター数をグラフから選別

>TS_log2.km<-kmeans(TS_log2[1:1000,1:6],14) #いよいよクラスタリング．引数の14はGap統計量で導出したクラスタ数

>TS_log2.km$cluster

>clusplot(TS_log2[1:1000,1:6],TS_log2.km$cluster,color = TRUE, shade = TRUE,labels=2,lines = 0) #結果のグラフ化。引数については参照先

>data.frame(TS_log2[1:1000,], TS_log2.km$cluster) #解析に用いた遺伝子行列にクラスター番号列を追加

>cluster_no<-data.frame(TS_log2[1:1000,], TS_log2.km$cluster) #上の行列を新しいオブジェクトに指定

>write.csv(cluster_no,"任意のパス") #csvファイルに書き出し、クラスターごとに発現変動折れ線グラフを書いてみる