チュートリアル初級では GC skew の解析を、既存のメソッドを利用することで簡単に行えることを説明しました。この中級編では実際にG-language GAEを用いて、GC skewの解析をする方法について説明します。

それでは早速プログラムを作ってみましょう。まずは前回もやったように、G-language をスクリプト内で起動します。

use G; my $gb = new G("ecoli"); G-language GAE で解析のためのプログラムを作る場合、ゲノムデータベースの読み込みやパース、そしてそれを操作する関数などについて時間を煩わされることは全くありません。上記の二行でおしまいです。つまり、解析自体に専念することが可能です。

そして今回GC skew のプログラムを作るわけですが、再利用可能なようにこれをサブルーチンとして作ることにしましょう。そして内部処理をしながら、最終的に他の解析でも使えるようにGC skewの値を返してくれるサブルーチンに設計しておくと便利です。第二引数として、GC skewの計算に必須なウインドウを与えるようにすれば汎用的になります。

  sub gcskew{ 
     my $gb = shift; 
     my $window = shift; 
     my @gcskew = (); 
 
     return @gcskew; 
  }

GC skew の計算は非常に簡単です。ゲノムを指定したウインドウに分け、それぞれのウインドウにおいて (C-G) / (C+G) を計算し、この推移をプロットするだけです。まずは外部のグラフ作成アプリケーションでこのプロット図がかけるように、カンマ区切りのテキスト形式である csv 形式で結果を出力するプログラムにしましょう。

G-language インスタンスでは、ゲノムの配列情報は $gb→{SEQ} に格納されています。このように、$gb→{LOCUS}, $gb→{BASE_COUNT} のような形で階層的にゲノムの情報が整理されているのです。

ウインドウ毎に、という処理には for 文が適しているでしょう。ゲノム全体の長さは length ($gb→{SEQ}), 文字のカウントは $a = $sequence =~ tr/a/a/、で出来ますね。

プログラムは以下のようになるはずです。

  my @location = ();  #Window 位置記憶のための配列 
  my $i = 0;  #カウンター初期化 
  for ($i = 0; $i * $window < length($gb->{SEQ}); $i ++){ # Windowがゲノムを超すまで繰り返す 
        my $sequence = substr($gb->{SEQ}, $i * $window, $window);  # Window内の配列を切り出す 
        my $c = $sequence =~ tr/c/c/;  #Cの数をカウント 
        my $g = $sequence =~ tr/g/g/;  #Gの数をカウント 
        my $skew = ($c-$g)/($c+$g);    #GC skew を計算 
        push (@location, $i * $window); #Window位置を保存 
        push (@gcskew, $skew);         #GC skewを保存 
  }

ここまで出来ればファイル出力は簡単ですね。

  my $j = 0; 
  open(OUT, '>gcskew.csv'); 
  print OUT "location,GC skew¥n"; 
  for ($j = 0; $j <= $i; $j++){ 
       print OUT $location[$j], ",", $gcskew[$j], "¥n"; 
  } 
  close(OUT);

さて、いとも簡単にGC skewを計算するプログラムは完成してしまいましたが、結果をファイルに出すのではそれを再加工しないと目で見える形での結果にはなりません。解析を効率化するのがG-language GAE ですから、自動的にグラフを作ってやりたいものです。

もちろんG-language GAEにはそんな時の為に簡単にグラフを作成するメソッドも備わっています。grapher()です。

  grapher(¥@x-axis, ¥@value1, ¥@value2 ...)
  option	description
  x	 x-axis label
  y	 y-axis label
  x1, x2, ...	 @value1, @value2, ... のラベル
  filename	 Filename
  title	 title

他にも豊富なオプションがありますが、詳しくはマニュアルを見てください。

つまりGC skew をグラフするだけならたったこれだけです。

grapher(¥@location, ¥@gcskew, -x⇒'bp', -y⇒'GC skew', -title⇒'GC skew', -filename⇒'gcskew.png'); G-language GAE 付属の簡易ビューア gimv で開くようにすれば解析と同時に結果が gnuplot で色鮮やかに描かれたものが表示されます。

  msg_gimv("graph/gcskew.png");

  use G; 
 
  my $gb = new G("ecoli"); 
 
  my @gcskew = &gcskew($gb, 10000); 
 
  sub gcskew{ 
     my $gb = shift; 
     my $window = shift; 
     my @gcskew = (); 
 
     my @location = ();  
     my $i = 0;   
     for ($i = 0; $i * $window < length($gb->{SEQ}); $i ++){      
           my $sequence = substr($gb->{SEQ}, $i * $window, $window);   
           my $c = $sequence =~ tr/c/c/;  
           my $g = $sequence =~ tr/g/g/;   
           my $skew = ($c-$g)/($c+$g);    
           push (@location, $i * $window);  
           push (@gcskew, $skew);      
     } 
 
     my $j = 0; 
     open(OUT, '>gcskew.csv'); 
     print OUT "location,GC skew¥n"; 
     for ($j = 0; $j <= $i; $j++){ 
           print OUT $location[$j], ",", $gcskew[$j], "¥n"; 
     } 
     close(OUT); 
 
 
     grapher(¥@location, ¥@gcskew, -x=>'bp', -y=>'GC skew', 
                              -title=>'GC skew', -filename=>'gcskew.png'); 
 
     msg_gimv("graph/gcskew.png"); 
     return @gcskew; 
  }