Tutorial-GCskew

Introduction

GC skew とは、一本鎖DNA分子におけるG含量とC含量のバイアスを表す指標で、(C-G)/(C+G) の式で表されます。Chargaffの塩基分配法則の第二項は一本鎖DNA分子においてG含量とC含量はほぼ等しくなると定めていますが、これはゲノム全体においての現象で、細かく区切られた領域ではバイアスが見られます。それどころか、バクテリアの種によっては綺麗にその傾向がシフトする個所がみられ、さらにはその個所が複製開始・終結地点と一致することが多いことが知られています。GC skew という現象が発生する原因には主に二説あり、それらはリーディング鎖とラギング鎖の異なる突然変異確率、そしてコドン使用による変異のバイアスなどによる、としています。

reference:
Lobry, J.R., 1996. Asymmetric substitution patterns in the two DNA strands of bacteria. Mol. Biol. Evol. 13, 660-665

Step 0 -G-language System の起動

それでは、GC skew の解析を実際に行う前にまず G-language System を起動しましょう。G-language System を使用すればゲノム解析は非常にシンプルです。例えばE.coli のGenbank形式のコンプリートゲノムを使用して解析をするなら、テキストエディタに以下の二行を書き込むだけで準備は終了します。

filename: test.pl

  use G; 
  $gb = new G("ecoli.gbk");

ためしにこのスクリプトを実行してみましょう。以下のような出力がでてくるはずです。

$ perl test.pl [ENTER] 
 
  Length of Sequence :   4639221 
           A Content :   1142136 (24.62%) 
           T Content :   1140877 (24.59%) 
           G Content :   1176775 (25.37%) 
           C Content :   1179433 (25.42%) 
              Others :         0 (0.00%) 
          AT Content :    49.21% 
          GC Content :    50.79%

この塩基含量の統計値の出力によって、G-language System は無事ファイルを読み込んだことを知らせています。

Genbank以外のゲノムデータベース（Fasta、EMBL、swiss、SCF、PIR、GCG、raw、aceなど）を使用したい場合は、

   use G; 
   $gb = new G("ecoli.fasta", "Fasta");

のように第二引数としてデータベースフォーマットを指定すれば読み込まれます。

Step 1 -GC skew の観察

G-language System にはゲノム解析のための多様な関数が標準で用意されています。幸い GC skew の解析もそんな標準関数の一つです。Step 0で作成したスクリプトにGC skew の関数を使う行を書き加えて実行してみましょう。

filename: test.pl

  use G; 
  $gb = new G("ecoli.gbk"); 
  gcskew(\$gb->{SEQ});

すると、以下のようなグラフが表示されるはずです。

これがE.coli ゲノムのGC skewを、10k bp のウインドウをずらして計算したものをグラフにしたものです。1,500,000 bp 周辺と3,800,000 bp 周辺にシフトポイントを持つ、非常に特徴的な傾向が見られることがわかります。

G-language System 起動時にロードされたゲノムデータは、全て $gb の中に格納されています。例えば全塩基配列は

  $gb->{SEQ}

の中にはいっています。G-language System の標準関数の多くはこの$gb、もしくは$gb->{SEQ}のリファレンス（\$gb->{SEQ} ：リファレンスがわからない人はワンセットで覚えてください）を引数として渡すことで動作します。gcskew()関数では$gb->{SEQ}のリファレンスが第一引数です。

Step 2 -GC skew の解析

G-language System の関数は容易に使えるだけでなく、多彩なオプションをもっています。

例えば gcskew()関数は以下のようなオプションを持ちます。

option description

-window GC skew を計算するウインドウサイズ（デフォルトは 10,000 bp)

-at 1 の時 AT skew を調べる（デフォルトは0）

-output f の時ファイル出力、g の時グラフ出力、show の時グラフ表示（デフォルトはshow）

-filename 出力ファイル名を指定。（デフォルトは-output がfの時 data/gcskew.csv、-outputがgの時graph/gcskew.gif）

-application 画像を表示するアプリケーション（デフォルトはgimp）

オプションは、

  gc(\$gb->{SEQ}, -window=>50000, -at=>1);

のように "-" をオプション名の頭に付け、"=>" で値と結びます。

それではより詳しくE.coli のGC skewを分析するために、E.coli のAT skew とGC skew を50k bp ごとのウインドウで見てみることにしましょう。AT skew とは GC skew 同様A含量とT含量のバイアスを見るもので、この傾向はGC skewと異なったり、GC skewほど顕著でなかったりします。

スクリプトは以下のようになります。

filename: test.pl

  use G; 
  $gb = new G("ecoli.gbk"); 
  gcskew(\$gb->{SEQ}, -window=>50000, -filename=>"gcskew50k.gif"); 
  gcskew(\$gb->{SEQ}, -window=>50000, -at=>1, -filename=>"atskew50k.gif");

グラフは以下の通りです。

GC skew: 50 k bp window

AT skew: 50 k bp window

50 k bp にするとGC skew のシフトポイントはより明確になりました。しかし、AT skew の結果はどうやらGC skewよりも顕著ではないようです。

課題2：いくつかの種で、GC skew とAT skew の傾向がどのように違うかを見てみましょう。また、ウインドウサイズを変えることで何か違ったものが見えてくるか試してみましょう。

Step 3 -GC skew を様々な角度から見る

Introductionで述べたように、GC skew の発生原因はリーディング鎖・ラギング鎖の突然変異の入り方の違い、そしてコドン使用によると考えられています。

リーディング鎖・ラギング鎖の突然変異の入り方の違いによってGC skew が発生したのなら、逆にGC skew から複製開始・終結地点を定義することが可能なはずです。実際GC skew のシフトポイントは頻繁に複製開始・終結地点と一致します。そこで、GC skew から複製開始・終結地点を予測してみます。

シフトポイントを特定するには累積GC skewが役に立ちます。累積GC skewはウインドウごとのGC skew値を累積するもので、シフトポイントが顕著に表れるのが特徴です。G-language Systemではこの累積GC skewは cum_gcskew() という関数になっています。

cum_gcskew()関数は以下のようなオプションを持ちます。

option description

-window Cumulative GC skew を計算するウインドウサイズ（デフォルトは 10,000 bp)

-at 1 の時 AT skew を調べる（デフォルトは0）

-output f の時ファイル出力、g の時グラフ出力、show の時グラフ表示（デフォルトはshow）

-filename 出力ファイル名を指定。（デフォルトは-output がfの時 data/cum_gcskew.csv、-outputがgの時graph/cum_gcskew.gif）

-application 画像を表示するアプリケーション（デフォルトはgimp）

   cum_gcskew(\$gb->{SEQ}, -window=>50000);

とすれば実行できますね。例えば以下のように表示されたはずです。

それではウインドウサイズを変えて累積GC skew のグラフをいくつか出力してみましょう。すると、シフトポイントが明確にあらわれることがわかるのではないでしょうか。

それでは累積GC skewを利用して複製開始・終結地点を予測してみましょう。G-language System には標準で find_ori_ter()という関数があり、この関数は内部で累積GC skew を使って複製開始・終結地点を予測します。

find_ori_ter()のオプションは　-window のみです。これは累積GC skewをとるウインドウの大きさで、デフォルトは1000になっています。この値が小さければ小さいほど予測値は精度が高くなりますが、逆に小さすぎると累積GC skew自体の計算が正確でなくなってしまい、結果的に精度が下がります。cum_gcskew()とあわせてウインドウサイズを変化させ、適切な複製開始・終結地点を予測してみましょう。

   find_ori_ter(\$gb->{SEQ}, -window=>500);

のように使います。

ウインドウサイズが10000bpなら、E.coliにおいて例えば以下のような予測が出力されます。

find_ori_ter: 
   Window size = 10000 
   Predicted Origin:   3915000 
   Predicted Terminus: 1545000

次に、GC skew の発生はコドン使用とかかわりがある、といわれていますから、GC skewがゲノム全体、遺伝子領域内、非コード領域、また、コドンの3塩基のうち最も「ゆらぎ」が大きいとされる三番目の塩基でどのようになっているか、比較して分析してみると面白い結果が見られるかもしれません。G-language System ではこの解析も genomicskew()という標準関数として備わっています。

genomicskew()関数は以下のようなオプションを持ちます。

option description

-divide GC skew を計算するウインドウの総数（デフォルトは 250個)

-at 1 の時 AT skew を調べる（デフォルトは0）

-output f の時ファイル出力、g の時グラフ出力、show の時グラフ表示（デフォルトはshow）

-filename 出力ファイル名を指定。（デフォルトは-output がfの時 data/genomicskew.csv、-outputがgの時graph/genomicskew.gif）

-application 画像を表示するアプリケーション（デフォルトはgimp）

   genomicskew($gb, -divide=>250);

として実行してみましょう。第一引数が　$gb であることに注意してください。

出力は以下のようになります。これまでの解析とあわせて考察してみましょう。

課題3：　GC skew の解析はまだまだ多角的に進めることが可能です。G-language System の多彩な標準関数を使用して、様々な解析を行ってみましょう。他種との比較などをしてみると興味深い結果がでるかもしれません。

Step 4 -より高度な解析のために

さて、G-language System には多数のゲノム解析関数が標準で用意されている為、それら関数のみでも幅広い解析を行うことが可能です。しかし、バイオインフォマティシャンとして生命現象の真理を探究するにはこれら関数をあくまで道具として利用し、より深く掘り下げて研究していく必要があります。

G-language System はゲノム解析の為の関数だけでなく、ゲノムデータベースを扱いやすくするためのプラットフォームを提供します。詳しくはG.pm のperldocドキュメンテーションに詳しく記されていますが、そのプラットフォームとは、$gb というG-language System　のインスタンスから呼び出せる関数であり、遺伝子ごとの処理やスタートコドン・ストップコドン周辺の処理、イントロン・エキソン処理、などです。

GC skew の解析も、それが最終目標ではなく途中で必要なツールになる場合が数多くあります。

例えば、遺伝子ごとに発現量とセカンドコドンの関係を調べ、その周辺の傾向をつかむためにGC skewを利用する、といった状況が想定されます。一例としてこれをどのようにG-language System を用いて表現するかを記します。

   use G; 
   $gb = new G("ecoli.gbk"); 
   cai($gb); 
   $w_val = w_value($gb); 
 
   foreach $cds ($gb->cds()){ 
      $secondcodon = $gb->after_startcodon($cds, 3); 
      $w_second = $$w_val{$secondcodon}; 
      $cai = $gb->{$cds}->{cai}; 
      if ($w_second < 0.5 && $cai > 0.8){ 
         $afterstart = $gb->after_startcodon($cds, 99); 
         gcskew(\$afterstart, -window=>9, -filename=>"$cds-gcskew.gif"); 
      } 
    }

まずG-language Systemを起動し、ゲノムデータベースをロードします。cai()関数を使い、G-language System　インスタンスにCAI値（Codon Adaptation Index：翻訳効率の指標だが遺伝子発現量の指標としても使われる)を付加します。コドン使用の偏りを示す値、W値をw_value()により取り出し、$w_valに格納します。

   foreach $cds ($gb->cds()){ 
 
   }

はG-language Systemを使いCDS領域毎の処理をするためのもっとも基礎的な方法です。$gb->cds()は$gbゲノムデータベースにある全てのCDSの名前を返します。つまり、それをforeachすれば全CDSについて解析することができるようになります。

$gb は

LOCUS, HEADER, FEATURE1 . FEATURE2 ... FEATURE4000 ... , CDS1 . CDS2 ... CDS4000 ..., SEQ

などの構造体を持ち、それぞれにFEATURE情報などが格納されています。つまり、各CDSの情報は　CDS＋番号　という名前の構造体に入っており、$gb->{CDS534}->{start} のようにそれぞれの情報に階層的にアクセスします。

      $secondcodon = $gb->after_startcodon($cds, 3); 
      $w_second = $$w_val{$secondcodon}; 
      $cai = $gb->{$cds}->{cai};

この部分は、まず標準関数であるafter_startcodon() で$cdsのスタートコドン以降3文字をとりだし、$secondcodonに入れています。また、そのW値を取得し、その遺伝子のCAI値も取得します。

      if ($w_second < 0.5 && $cai > 0.8){ 
         $afterstart = $gb->after_startcodon($cds, 99); 
         gcskew(\$afterstart, -window=>9, -filename=>"$cds-gcskew.gif"); 
      }

この部分では、CAI値が0．8以上、つまり発現量が非常に高い遺伝子でスタートコドン直後にW値が0.5未満のコドン、つまりレアコドンが使用されている場合、スタートコドン下流99bp のGC skew をウインドウ9bp (つまり3コドン) ずつ見る、というものです。出力オプションはデフォルトで「表示」なので、つまり該当する遺伝子があるだけグラフが自動的に表示されることになります。これらの出力をもとに、特異な発現機構を発見することが可能かもしれません。

ここに示した例はある程度複雑ですが、それでも実態はわずか13行のスクリプトです。G-language Systemではこのように複雑な解析を効率よく、容易に行うことを可能にします。

それでは、あとはこのG-language Systemを使ってあなたが生命の真理を探究してみてください。

option	description
-window	GC skew を計算するウインドウサイズ（デフォルトは 10,000 bp)
-at	1 の時 AT skew を調べる（デフォルトは0）
-output	f の時ファイル出力、g の時グラフ出力、show の時グラフ表示（デフォルトはshow）
-filename	出力ファイル名を指定。（デフォルトは-output がfの時 data/gcskew.csv、-outputがgの時graph/gcskew.gif）
-application	画像を表示するアプリケーション（デフォルトはgimp）

option	description
-divide	GC skew を計算するウインドウの総数（デフォルトは 250個)
-at	1 の時 AT skew を調べる（デフォルトは0）
-output	f の時ファイル出力、g の時グラフ出力、show の時グラフ表示（デフォルトはshow）
-filename	出力ファイル名を指定。（デフォルトは-output がfの時 data/genomicskew.csv、-outputがgの時graph/genomicskew.gif）
-application	画像を表示するアプリケーション（デフォルトはgimp）