スライド 1

Size: px

Start display at page:

Download "スライド 1"

もえりかせ
4 years ago
Views:

1 遺伝子発現データの階層的非階層的クラスタリング慶應義塾大学環境情報学部

2 階層的クラスタリング (1) Data1 Data2 Data3 Data 相関行列 Data1 Data2 Data3 Data4 Data Data Data Data Data1 Data2 Data3 Data4

3 階層的クラスタリング (2) Data1 Data2 Data Data1 Data2 Data3+4 Data Data Data Data1 Data2 Data3 Data4

4 階層的クラスタリング (3) Data1+2 Data Data1+2 Data3+4 Data Data Data1 Data2 Data3 Data4

5 ヒト繊維芽細胞の発現データのクラスタリング細胞修復血管形成細胞周期コレステロールの合成 time 早期応答 Eisen M et al. 1998

6 階層的クラスタリングのステップ 1. 全クラスター間の相関係数の計算 2. 最も相関が高いクラスターを併合 3. 複数のクラスターがあるなら 1. に戻る

7 クラスター間の相関係数はどのように計算する? クラスター A, Bの要素 a, bの相関係数をr a,b として最短距離法 (Single Linkage Method) max(r a,b ) 最長距離法 (Complete Linkage Method) min(r a,b ) 群平均法 (Average Linkage Method) average(r a,b )

8 階層的クラスタリング (Perl) # クラスター a, b の要素間の相関を計算する sub corr_clust($$$){ my $clust_a = shift; my $clust_b = shift; my $d = shift; my $max = -1.0; for my $i (@$clust_a){ for my $j (@$clust_b){ if($max < $d->[ $i ]->[ $j ]){ $max = $d->[ $i ]->[ $j ]; return $max; # 次元数 dim のデータ群 data を階層的クラスタリング sub h_cluster($){ my $data = shift; # データセット $data->[ データ番号 ]->[ データの次元番号 ] # クラスター情報 $cl[ クラスター番号 ] = [ データ番号,,, ] # クラスター情報の初期化 foreach(0..$#$data){ $cl[$_] = [ $_ ]; # 相関行列の作成 for my $i (0..$#$data){ for my $j (0..$#$data){ $d[$i]->[$j] = corr_single($data->[$i], $data->[$j]); # 階層的クラスタリング while ($#cl > 0){ my $corr_max = -1.0; my($c1, $c2); for (my $i = 0; $i < $#cl; $i ++){ for (my $j = $i + 1; $j <= $#cl; $j ++){ my $corr = corr_clust($cl[$i], if ($corr > $corr_max){ $corr_max = $corr; $c1 = $i; $c2 = $j; printf("[ Merge %d ] Merging the following two clusters: n", $#$data - $#cl + 1); printf("(1) Cluster containing data #"); print " n"; printf("(2) Cluster containing data #"); print " n"; printf("distance between these two clusters = %lf n n", $corr_max); $cl[ $c1 ] = merge($cl[ $c1 ], $cl[ $c2 ]); # クラスター c2 を c1 へ統合して格納 splice(@cl, $c2, 1) # クラスター c2 を破棄クラスターの統合クラスター間の相関係数

9 階層的クラスタリング (Perl) その他の設定 use List::Util qw(sum); # n 次元データ a, b の相関係数を計算する sub corr_single($$){ my $a = shift; my $b = shift; my $n = $#$a + 1; my $mean_a = sum(@$a) / ($#$a + 1); my $mean_b = sum(@$b) / ($#$b + 1); my $var_a = 0; my $var_b = 0; foreach(0..$n-1){ $var_a += ($a->[$_] - $mean_a)**2 / $n; foreach(0..$n-1){ $var_b += ($b->[$_] - $mean_b)**2 / $n; my $corr = 0; foreach(0..$n - 1){ $corr += (($a->[$_] - $mean_a) / ($var_a ** 0.5)) * (($b->[$_] - $mean_b) / ($var_b ** 0.5)) / $n; return $corr;

10 階層的クラスタリング (C 言語 ) /* 次元数 dim のデータ群 data を階層的クラスタリングデータ数 N_DATA はグローバル変数として与えられる */ void h_cluster(double data[n_data][dim]){ int i, j, n; double d[n_data][n_data]; /* データ間の相関行列 */ double corr, corr_max; /* データ間の相関 */ int c1, c2; static struct CLUSTER cl[n_data]; /* クラスター情報 */ static struct CLUSTER w; /* ワーク用 */ /* クラスター情報の初期化 */ for (i=0; i < N_DATA; i++){ cl[i].n = 1; cl[i].gene[0] = i; /* 相関行列の作成 */ for (i = 0; i < N_DATA; i++) for (j = 0; j < N_DATA; j++) d[i][j] = corr_single(data[i], data[j], DIM); /* 階層的クラスタリング */ n = N_DATA - 1; while (n > 0){ corr_max = -1.0; for (i = 0; i < n; i ++){ for (j = i + 1; j <= n; j ++){ corr = corr_clust(cl[i], cl[j], d); if (corr > corr_max){ corr_max = corr; c1 = i, c2 = j; printf("[ Merge %d ] Merging the following two clusters: n", N_DATA - n - 1); printf("(1) Cluster containing data #"); for(i = 0;i < cl[ c1 ].n;i ++){ printf("%d", cl[ c1 ].gene[i]); if(i < cl[ c1 ].n - 1)putchar(','); printf(" n(2) Cluster containing data #"); for(i = 0;i < cl[ c2 ].n;i ++){ printf("%d", cl[ c2 ].gene[i]); if(i < cl[ c2 ].n - 1)putchar(','); printf(" ndistance between these two clusters = %lf n n", corr_max); merge(cl[ c1 ], cl[ c2 ], w); /* クラスターを統合 */ cl[ c1 ] = w; /* クラスター c2 を c1 へ統合して格納 */ cl[ c2 ] = cl[n]; /* クラスター c2 を破棄し一番後ろのクラスターを代わりに格納 */ n--; /* クラスターを 1 つ減らす */

11 階層的クラスタリング (C 言語 ) その他の設定定数の定義構造体の宣言など #define N_DATA 4 /* データ数 */ #define DIM 21 /* クラスター情報 */ struct CLUSTER { int n; /* データ数 */ int gene[n_data]; /* 保持しているデータ番号の集合 */ ; 相関係数の計算 /* n 次元データ a, b の相関係数を計算する */ double corr_single(double a[], double b[], int n){ int i; double d = 0.0; double mean_a, mean_b; double var_a, var_b; double corr; for(mean_a = 0, i = 0;i < n;i ++)mean_a += a[i] / n; for(mean_b = 0, i = 0;i < n;i ++)mean_b += b[i] / n; for(var_a = 0, i = 0;i < n;i ++) var_a += (a[i] - mean_a)*(a[i] - mean_a) / n; for(var_b = 0, i = 0;i < n;i ++) var_b += (b[i] - mean_b)*(b[i] - mean_b) / n; for(corr = 0, i = 0;i < n;i ++) corr += (a[i] - mean_a) / sqrt(var_a) * (b[i] - mean_b) / sqrt(var_b) / n; return corr;

12 クラスタリングの計算時間対象がN 個あるなら以下の数の組み合わせの中から最も関連が高いものを選ぶ 1 N( N 1) 2 クラスターの併合を N=1 になるまで繰り返すなら調べる組み合わせの数は N i= i( i 1) = N( N 結局これは O(N 3 ) の計算量となりデータ数が多いと膨大な計算時間になる 2)

13 非階層的クラスタリング K 平均法自己組織化マップ

14 発現データを二次元平面へ発現レベル 5 A 条件 B B A 条件 1 条件条件 1

15 K 平均法 1. 参照点の設定 2. 各要素の所属を最も近い参照点とする 3. 各クラスターの重心を計算 4. 参照点を重心の方向へわずかに移動させるに戻る

16 所属参照点の決定 (a) 5 10 (b)

17 重心の決定と参照点の移動 (c) 5 10 (d)

18 繰り返し (e) 5 10 (f)

19 K 平均法デモ

20 Perl による K 平均法の実装の本日の授業の箇所よりファイルをダウンロードし解凍ファイルを各自の作業ディレクトリに配置 kmeans_func.pm k 平均法を実装する上で有用な関数を集めたファイル Plot_field.pm, Plot_field_KM.pm 2 次元データを画面上にプロットするときに有用な関数を集めたファイル kmeans_test.pl 上記ファイル中に含まれる関数の使用法のサンプル testfile.txt 遺伝子発現の 4 次元架空データ testfile2d.txt 遺伝子発現の 2 次元架空データ eisen_exp.txt 酵母菌の遺伝子発現の実データ

21 初期設定 kmeans_func.pm を同一ディレクトリに配置プログラムを作成する但し先頭には以下のように記述 #!/usr/bin/perl -w use strict; require Exporter; use kmeans_func qw(determine_cluster calc_gravity_center initialize); *::CLUSTERS = 3; # クラスター数

22 要素の数次元数データの読み込み遺伝子名機能 my($points, $DIM, $gene) $ARGV[0]); クラスター数入力ファイル参照点の位置要素の位置遺伝子名入力ファイルの中味の例 (TAB 区切り ) 遺伝子の機能 $gene->{$i->{ name が $i 番目の遺伝子の名前 $gene->{$i->{ func が $i 番目の遺伝子の機能 $points[$i]->[$j] で $i 番目の要素の $j 次元目の値を表す $ref_position[$i]->[$j] で $i 番目の参照点の $j 次元目の値を表す発現データ Gene A Transcription Gene B Translation Gene C Translation Gene D Transport Gene E Transport Gene F Transport Gene G Transcription

23 my($points, $DIM, $gene) $ARGV[0]); print "***** Read positions from file $ARGV[0] ***** n n"; 実装例 for(my $i = 0;$i < $POINTS;$i ++){ print "Point #$i t", $gene->{$i->{ "name"," - "; print $gene->{$i->{ "func"," t"; print join(" " n"; print " n"; print "***** Generated reference positions ***** n n"; for(my $i = 0;$i < $::CLUSTERS;$i ++){ print "Cluster #$i t"; for(my $j = 0;$j < $DIM;$j ++){ printf("%+.3lf ", $ref_position[$i]->[$j]); print " n"; 実行結果 ***** Read positions from file testfile.txt ***** Point #0 Gene A - Transcription Point #1 Gene B - Translation Point #2 Gene C - Translation Point #3 Gene D - Transport Point #4 Gene E - Transport Point #5 Gene F - Transport Point #6 Gene G - Transcription ***** Generated reference positions ***** Cluster # Cluster # Cluster #

24 所属クラスターの決定各要素の所属各要素の参照点の座標 @points); クラスター $i に要素 $j が所属するなら $points_belong[$i]->[$j] = 1 そうでなければ $points_belong[$i]->[$j] = 0

25 実装例 print "***** Belongings of each points ***** n n"; for(my $i = 0;$i < $::CLUSTERS;$i ++){ print "Cluster #$i: "; for(my $j = 0;$j < $POINTS;$j ++){ printf("%d ", $points_belong[$i]->[$j]); print " n"; print " n"; 実行結果 ***** Belongings of each points ***** Cluster #0: Cluster #1: Cluster #2: 要素 0 要素 2 要素 1

26 重心の計算要素の所属クラスター各要素の参照点の座標各クラスターの重心

27 @target_position, 実装例 print "***** Gravity center positions ***** n n"; for(my $i = 0;$i < $::CLUSTERS;$i ++){ print "Cluster #$i: "; for(my $j = 0;$j < $DIM;$j ++){ if(defined($target_position[$i]->[$j])){ printf("%+.3lf ", $target_position[$i]->[$j]); else { print " "; print " n"; print " n"; 実行結果 ***** Gravity center positions ***** Cluster #0: Cluster #1: Cluster #2:

28 結果を端末で分かりやすく図示 ( 二次元データのみ ) Plot_field.pm および Plot_field_KM.pm を同一ディレクトリに配置 use Plot_field_KM; my $fld = new Plot_field_KM -10, -10, 10, $::CLUSTERS);

29 (-10,10)...(10,10) >0< >2< [0]......[2] [1] >1< (-10,-10)...(10,-10) 各要素とその所属クラスター ( この例ではクラスター 0に所属 ) クラスター 0 の重心クラスター 0 の参照点

30 C++ による K 平均法の実装の本日の授業の箇所よりファイルをダウンロードし解凍ファイルを各自の作業ディレクトリに配置 Makefile コンパイル方法が書かれているファイル kmeans_mac.c プログラムの main このファイルを改良して K 平均法を実装する point.? 各要素のクラスおよびそれに対する有用な関数を集めたファイル plot_field.? クラスタリングの様子をグラフィカルに描画するためのファイル rsgraph.? X11 を操作する関数を集めたファイル initialize.? 各要素や参照点の初期位置などをユーザに求める関数以下のように実行 make./kmeans

31 K 平均アルゴリズムのメインループ for(int k = 0;k < 100;k ++){ determine_cluster(points_belong, target_position, ref_position, points, n_points, n_refs); /* 所属クラスターの決定 */ calc_gravity_center(points_belong, target_position, ref_position, points, n_points, n_refs); /* 重心の計算 */ fld.plot_kmeans(points_belong, target_position, ref_position, points, n_points, n_refs, k); /* 描画 */ /* Here, write your own code */

32 メインループの各変数の意味 class point に属するオブジェクトは四則演算が可能 target_position は現在のクラスターの重心 ref_position は参照点の座標 points が各要素の座標 n_points は要素数 n_refs は参照点の数 points_belong[i][j] = 1 なら要素 i はクラスター j に属している属していないなら 0

33 プロモータ配列の探索発現パターン TF A B C 共通の配列

34 発現データのクラスタリングとプロモータ探索

35 プロモータ配列の探索発現データのクラスタリング顕著な上流配列パターン顕著な上流配列パターンを持つクラスター ( 配列パターンの特異性の解析 )

36 発現パターンが類似している遺伝子の翻訳産物 ( タンパク質 ) は相互作用する? TF A 発現パターン B 発現パターンが類似

37 相互作用するペアとしないペアの遺伝子発現の相関係数の分布 (Grigoriev. 2001)

38 K 平均法による解析 (Ge et al. 2001)

39 K 平均法による解析 (Ge et al. 2001)

40 新規複合体サブユニット (Ge et al. 2001)

スライド 1

スライド 1 遺伝子発現データの階層的クラスタリング慶應義塾大学環境情報学部遺伝子の発現のプロセス DNA ATG 遺伝子 TAA 転写 mrna AUG UAA 翻訳タンパク質遺伝子はいつ発現している? 全ての遺伝子が常に発現しているわけではない皮膚の細胞と肝臓の細胞の DNA 情報は同一細胞の形質が異なるのは発現している遺伝子が異なるため ON OFF 細胞の分化 (1) 細胞が不可逆的に特殊化していくこと