ex603 - nibb-gitc/gitc2025mar-rnaseq GitHub Wiki

ex603

(発展問題: 教師ありクラスタリング手法であるk-means法はトレーニングコースでは解説していません)

事前にクラスター数を指定して行うことを「教師あり」クラスタリングと呼ぶ。k-means法は教師ありクラスタリングの1種で、汎用される手法である。 教師ありクラスタリングは前提(仮説)が正しい時には教師なしクラスタリングよりも精度よくクラスターを特定できる。 また、この例には含めていないが、k-means法ではクラスターの重心プロファイル(イメージとしては計算したクラスターの平均的な遺伝子発現プロファイル)を計算できる他、各サンプルを重心プロファイルからどの位、離れているか(似ていないか)という特徴量エンジニアリングにも使用できる。

この演習問題で使用する変異体(正確には変異体と野生型のトランスクリプトームのlog2 fold change)が関与するシグナル伝達経路は

  • coi1, dde2, jar1, jin1: ジャスモン酸シグナル伝達経路
  • ein2-1, ein3: エチレンシグナル伝達経路
  • npr1-1, sid2-2: サリチル酸シグナル伝達経路

である。これらの変異体遺伝子発現プロファイルをk-means法を使って解析し、クラスター数の妥当性を考察せよ。kの数を1から7に変化させて解析し、クラスター数を評価しなさい。

inputMatrix <- read.delim(
  "gitc/data/MS/Sato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_small.txt", 
  header=TRUE, 
  row.name=1
)
ET_JA_SA_input_matrix <- inputMatrix[,c("coi1", "dde2", "jar1", "jin1", "ein2.1", "ein3", "npr1.1", "sid2.2")]
kmeans_results        <- sapply(1:7, function(x) kmeans(t(ET_JA_SA_input_matrix), centers=x, iter.max = 100)$cluster)
print(kmeans_results)

[注意] centers引数、iter.max引数を調整することにより、結果が変化します。

       [,1] [,2] [,3] [,4] [,5] [,6] [,7]
coi1      1    2    3    2    1    1    2
dde2      1    1    2    4    5    5    5
jar1      1    1    2    3    5    3    4
jin1      1    1    2    3    5    3    3
ein2.1    1    1    2    3    3    4    6
ein3      1    1    2    3    5    3    4
npr1.1    1    1    1    1    4    6    1
sid2.2    1    1    1    1    2    2    7

k-means法はセントロイドの初期値がランダムに指定されるので、毎回同じ答えを返すとは限りません。よって、解析を繰り返して実行し、結果が収束しているかを確認する必要があります。

以下のスクリプトで解析結果の多数決を確認してみましょう。この例ではk=4で解析してみます。

make_k_means_clusters_to_string <- function(x, k){
  k_means_result <- kmeans(x, centers=k, iter.max = 100)$cluster
  sample_names   <- names(k_means_result)
  group_strings  <- sapply(1:k, function(x) 
                           paste(sort(sample_names[which(k_means_result == x)]), collapse=",")
                    )  
  return( paste(sort(group_strings), collapse=" : "))
}

k_means_repeat_results <- sapply(1:1000, function(x) make_k_means_clusters_to_string(t(ET_JA_SA_input_matrix), 4))
table(k_means_repeat_results)

今回は以下の結果になりました。

k_means_repeat_results
coi1 : dde2 : ein2.1,ein3,jar1,jin1 : npr1.1,sid2.2 coi1 : dde2,ein2.1,ein3,jar1,jin1 : npr1.1 : sid2.2 
                                                208                                                 215 
coi1 : dde2,ein3,jar1,jin1 : ein2.1 : npr1.1,sid2.2 
                                                577