User Tools

Site Tools


introductiontog-languagejapanese

1.はじめに

初めて一つの生物の完全ゲノム配列が決定された1994年からわずか10年あまりで,ゲノムの配列決定技術は飛躍的な進歩を遂げ,今ではヒトを含む数百種もの生物のゲノム情報が誰でもインターネットで入手できるようになっています. National Center for Biotechnology Information (NCBI) のGenBankデータベースに登録されている塩基配列情報はこれまでに指数関数的な増加を見せており,シーケンス技術の目覚ましい発展を見る限り,今後もこの勢いはとどまることがなさそうです.ゲノム配列決定とともに,網羅的データはトランスクリプトーム,プロテオーム,メタボローム,インタラクトームなど多岐に渡る階層で急速に蓄積されてきていますが,これらの定量データは生物学につきものである実験誤差や測定ノイズをある程度は含んでしまっており,擬陽性データを省き,扱うデータの信頼性を保つために,統計的なデータの前処理が必要になっています.また,これらのデータは網羅的とはいっても,完全に網羅していることは稀です.当然これらのデータが重要であることは間違いないですが,ゲノム情報は全遺伝情報のセットであり,本当の意味で網羅的である点,そして,いわばデジタルデータと言ってよいほどノイズが少ない点など,他のデータとは性質上一線を画しています.このようなゲノム情報特有の性質のため,今後も生物学においてゲノム情報は常に基軸として重要な位置を占め続けるでしょう.

一方で,増え続ける大量の配列情報を有効活用するためには,優れたコンピュータ解析環境が不可欠となります.数百種のバクテリアゲノムの比較解析を行うためには計算機の速度とともに効率の良い計算アルゴリズムも必要になり,同じような解析を多種にくりかえし適用する場合には,現存する様々な手法や道具を効果的に再利用できるような仕組みが必趁です.また,バイオインフォマティクスにおける研究は,生物実験におけるプロトコルのように,一連の手続きをワークフローとして繋げることによってはじめて成し得るものなので,一連の解析をパイプラインとして繋げてワークフローを構築できる環境も必要になります.本章では,このような目的を達成するソフトウェア環境であるG-language Genome Analysis Environmentを使い,特にこのソフトウェア環境が得意とするバクテリアゲノムの解析に関して,Perl言語によるプログラミングを交えながら具体的解析例を紹介していきます.

1.1 G-language GAE とは

G-language Genome Analysis Environment (G-language GAE)は,慶應義塾大学先端生命科学研究所で著者の一人である荒川が中心となって2001年から開発している,バイオインフォマティクスのための統合パッケージです.現在200以上の解析プログラムが含まれており,グラフィカルなユーザインタフェースからマウスの操作のみでもある程度の解析までは用意に行うことができます.より高度な解析のためには独自のコマンドラインの対話型シェルを備え,さらにプログラミング言語Perlのライブラリとして扱う事ができ,非常に強力な拡張性を持っています.そのため,Perl言語で書かれたその他のプログラムや,Comprehensive Perl Archive Network (CPAN)に登録されている豊富なライブラリやパッケージ,そしてBioPerlとシームレスに連動させることができます.このシステムはGNU General Public Licenseによって公開されているオープンソースのフリーソフトウェアなので,いつでもユーザはシステムのソースコードを見て,必要に応じて改変することができます.

200以上の解析プログラムは特にバクテリアゲノム解析に関するものが充実していますが,G-language GAEはゲノム情報解析にとどまらず,マイクロアレーやプロテオームの解析,さらにはそれらの情報を統合的に扱うシステムバイオロジー研究までを視野に入れた幅広い分野の解析を含んでいます.特に同じ研究所で開発されている細胞シミュレーションソフトウェアであるE-Cellとも連携できるため,細胞シミュレーションのためのデータマイニング環境として使う事もできます.このような解析に関しては本書では扱いませんので,プロジェクトのウェブページをご覧ください.

参考文献:

  • Arakawa. et al., 2003. G-language Genome Analysis Environment: a workbench for

nucleotide sequence data mining. Bioinformatics. 19, 305-306

  • Arakawa. et al., 2006. G-language System as a platform for large-scale analysis

of high-throughput omics data. Journal of Pesticide Science. 31, 282-288

2. G-language GAEの基本的な使い方

2.1 グラフィカルユーザインタフェースによる解析

G-language GAEには,基本的にはマウスのクリックのみで高度な解析を行うことができるグラフィカルユーザインタフェース(GUI)があります.GUIを起動するには,KNOBのメニュー(サイエンス&数学→G-language→glang)から選択するか、コマンドラインでglangとタイプします.

 knoppix@Knoppix:~$ glang

すると,3つのウインドウ:解析の実行などを行うメインのウインドウであるG-language Manager,システムメッセージを表示するConsole,そして解析の詳細を調節するConfigurator,が起動します.起動した状態で,すでにG-language GAEにはバクテリア解析システムというソフトウェアと,解析対象としてのマイコプラズマ菌のゲノム,そして解析を行うためのパラメータなどがセットアップされています.まずはそのままG-language ManagerのExecuteボタンを押して,初期設定のままで解析を行ってみましょう.

                      codon_usage

               base_information_content                       genomicskew

                      Text Output                                                view_cds

                                                           genome_map

実行すると,すぐに次々と結果が出力されてきます.文字列として出力される結果はText Outputウインドウに,それ以外の結果は個別のウインドウで表示されます.G-language GAEでは結果はできる限りグラフや表などにし,一目で結果を理解できるようにしているため,ほとんどの解析はグラフィカルな出力をします.これによりユーザは煩わしいデータの整理を最小限にとどめて研究や考察を進㒁ることができます.

G-language GAE全体では200以上の解析が実装されていますが,ここで使っているバクテリア解析システムは,特にバクテリアのコンプリートゲノムの解析に適している約30個の解析が用意されています.例えば,今回出力されたcodon_usageはコドンテーブルを表示するもの,base_information_contentはスタートコドン周辺領域の情報量を算出するもの,genomicsskewは様々な領域におけるGC skewを表示するもの,view_cdsはスタートコドン及びストップコドン前後100bpの塩基含量をグラフにしたもの,そして,genome_mapはゲノム中の塩基含量や遺伝子などの情報を地図として表示するものです.

2.1.1 G-language Manager

GUIのメインウインドウであるG-language Managerには,3つのボタンが用意されています.ConfigureボタンはConfiguratorウインドウを開き,Executeボタンは既に例で見たように,選択されている解析を実行します.Generateボタンは,選択されている解析を実行する代わりに,その内容をPerlのスクリプトとして出力します.出力されたPerlスクリプトを編集することで,より細かい設定を行ったり機能拡張を行うことができます.

Helpメニューでは,G-language GAEに関する情報を得たり,プロジェクトのウェブページを見ることができます.プロジェクトのウェブページには様々なドキュメントがある他,最新版のダウンロードなどができます.Fileメニューでは,アプリケーションの終了や,解析システムの読み込みができます.File→Load ScriptメニューではG-language Configuration File (GCF) 形式の解析システムファイルを選択して読み込むことができます.標準ではバクテリア解析システムのdefault.gcfが選択されていますが,ここでインストールパッケージのshare/gcfにある様々なGCFファイルを選ぶことによって幅広い解析を行うことができます.以下に標準で提供されている GCFファイルのリストを示します.プロジェクトのウェブページからは追加のGCFファイルを取得することができます.

GCFファイル名 説明
all.gcf G-language GAEに含まれる200以上の全解析プログラムを含む.
default.gcf バクテリア解析システム.バクテリアのコンプリートゲノムの解析に適した解析.
oasys.gcf 特定の遺伝子に関する情報を網羅的に計算するシステム.
CASYS.gcf cDNA解析システム.cDNAデータをクラスタリングし,選択的スプライシング候補をあげる.
Mapping.gcf cDNAやESTデータをゲノム上にマップするシステム.
STeP.gcf STS-ePCRシステム.ePCRを用いてOMIMデータなどをSTS間の場所にマップする.

File→Convert Perl Scriptメニューでは,Perlのプログラムを読み込んですぐにGUIに変換することができます.Perlプログラム内のサブルーチンの入出力を自動的にシステムが判別し,それをGCF形式に変換してGUIに読み込みます.File→Conovert Perl ScriptメニューとGenerateボタンを使うことにより,PerlスクリプトからGUIへ,そしてまたPerlスクリプトへ,と,インタフェースを繰り返し変えながらソフトウェアの改良を行うことが可能です.

2.1.2 Configurator

Configuratorでは解析に使うプログラムを選択し,それらの詳細なオプションを設定することができます.最上部の$gb部分に使用するデータベースのファイル名を指定し,左側にある解析プログラムのリストから,実際に使用するものを選びOn/Offボタンを押して使用/不使用を切り替えて行きます.On/Offボタンの右側の空白に数字を入力することで,解析を行う順序を定義することもできます.

より詳細な解析の設定を行いたい場合,任意の解析のEditボタンをクリックすると,Configuratorの右側部分にその解析プログラムの利用可能なオプション一覧が現れます.ここで変更したいオプションの値を入力することができます.一連の操作が終了し,使用するプログラムとそのオプションを設定し終えたら,Configurator下部のSaveボタンを押すと変更した設定が保存されます.あとはG-language ManagerのExecuteボタンを押せば解析がスタートします.

オプションは解析プログラムによって様々なものがありますが,多くに共通しているものとして,ゲノムや塩基配列などを示すinstance_G(通常は$gbそのままで変更の必要なし),ファイルに出力する場合のファイル名である-filename,返り値を格納する変数を指定する-Return,そして結果の出力形式を指定する-output(fはファイルへ出力,gはグラフへ出力,showはグラフへ出力後表示,stdoutは標準出力),があります.

2.1.3 解析をカスタマイズしてみましょう

Configuratorで,ゲノム中のGC含量の推移をプロットするgcwinのみをOnにして他を全てOffにしてみてください.その状態をSaveしてExecuteします.

次に,Editボタンを押してオプションを表示させ,windowを1000や5000などに変えてSave/Executeをしてみましょう. windowオプションはゲノムを区切るウインドウのサイズを指定するもので,デフォルトでは10000になっています.マイコプラズマ菌のゲノムは非常に小さいので,デフォルト値では大きすぎるのがわかるのではないでしょうか.

最適なウインドウがみつかったら,今度はatオプションを1にしてみましょう.このオプションを1にすると,GC含量ではなくAT含量の推移を見ることができます.

同様に,興味のある解析をOnにして,オプションをいろいろとカスタマイズして試してみてください.

 コラム: g2s とは
 G-language GAEのGUIを起動する時に,KNOBのメニュー(サイエンス&数学→G-language)にg2sというアプリ
 ケーションが見えたと思います.これは genome to shootingの略で,汎用ゲノム解析ソフトウェアパッケージG-
 languageの高度なゲノム解析能力をフルに使い,ゲノムの情報を元にシューティングゲームを自動生成するアプリケー
 ションです.
 
 時間の代わりにゲノムの塩基を進んでゆき,遺伝子領域に到達すると敵が出現!出現方向は二重らせんの鎖の方向,攻撃
 力はコドンバイアスから算出した発現量,種類はCOGオーソログデータベースにより検索・取得した遺伝子機能.ゲーム
 を遊びながら,複製方向による遺伝子方向バイアスや特定機能遺伝子のクラスター,さらには遺伝子名まで気づかぬうち
 に学習できます.エネルギーをためて放つ無敵の波動砲,回復系と攻撃系のアイテム,手に汗握るターボモードなど,
 ゲーム性も抜群!ゲノムによってステージが変わりますので,すでに数百あり,毎月数個以上新たに読まれるバクテリア
 ゲノムを使えば,ほぼ無限に楽しむことができます.研究の息抜きに是非遊んでみてください.
 
 操作方法は,十字キーで移動,スペースキーでミサイル発射(最大連続5発まで),シフトキーで全てを貫通する波動法
 発射(右のゲージで確認できるエネルギーを消費します),pキーでゲームの中断・再開ができます.敵を倒すとたまに
 落とす赤のカプセルで波動エネルギー充填,青いカプセルでライフ回復ができます.

2.2 G-language シェル

2.2.1 まずは起動して使ってみましょう

G-language GAEは,対話型で解析を進めることができるシェルを備えています.このシェル上ではPerlの全ての機能とG-language GAE関連の全ての機能,そして簡単なUNIXシェルの機能を使うことができ,非常に強力なインタフェースです.起動するにはGとコマンドラインで入力します.

   knoppix@Knoppix:~$ G
   
                __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/
                   
                      G-language  Genome Analysis Environment v.1.7.2
   
    
                                http://www.g-language.org/
   
                 Please cite: 
                    Arakawa K. et al. (2003) Bioinformatics.
                    Arakawa K. et al. (2006) Journal of Pestice Science.
   
                 License: GNU General Public License
                 Copyright (C) 2001-2007 G-language Project
                 Institute for Advanced Biosciences, Keio University, JAPAN 
   
                __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/
              
   G >

行頭に G > と表示されたら,続けてPerlによるプログラムや,サポートされているシェルコマンド(ls, cd, cp, rm, less, convert)を入力すればすぐに結果が表示されます.ファイル名に関してはPerlコマンド中であっても,タブによる補完機能があります.また,通常のUNIXシェル同様,上キーを押すことによってhistory機能により直前に入力したコマンドを呼び出すことができる他,Ctrl-Kで一行消去したりすることができます.

   
   G >$value = 256;
   
   G >print $value * 4;
   1024
   G >ls
   Desktop	  book   tmp
   G >$mgen = new G("mgen");
   
   
   Accession Number: NC_000908
   
     Length of Sequence :    580074
              A Content :    200543 (34.57%)
              T Content :    195695 (33.74%)
              G Content :     92312 (15.91%)
              C Content :     91524 (15.78%)
                 Others :         0 (0.00%)
             AT Content :    68.31%
             GC Content :    31.69%
   
   G >
   

シェルを終了するにはquitコマンドを使います.すると,ワークスペースを保存するかと聞いてきます.ここでyを入力すると,作業中に作成した全ての変数を永続的に保存しておいて,次回シェルを起動した時にすぐに利用可能なように呼び出してくれます.通常プログラムではメモリの内容は終了時に消去されますが,シェルの永続メモリ機能はプログラム終了後も維持され,次回起動時に自動的に読み込まれます.

   G >quit
     save value ?  y/n ?
   quit >>y
   knoppix@Knoppix:~$ G
   
                __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/
                   
                      G-language  Genome Analysis Environment v.1.7.2
   
    
                                http://www.g-language.org/
   
                 Please cite: 
                    Arakawa K. et al. (2003) Bioinformatics.
                    Arakawa K. et al. (2006) Journal of Pestice Science.
   
                 License: GNU General Public License
                 Copyright (C) 2001-2007 G-language Project
                 Institute for Advanced Biosciences, Keio University, JAPAN 
   
                __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/
              
        G >$mgen->seq_info()
        
        
   Accession Number: NC_000908
   
     Length of Sequence :    580074
              A Content :    200543 (34.57%)
              T Content :    195695 (33.74%)
              G Content :     92312 (15.91%)
              C Content :     91524 (15.78%)
                 Others :         0 (0.00%)
             AT Content :    68.31%
             GC Content :    31.69%
   
   
   G >
   

また,シェルで作業をした後にmakelogコマンドを入力すると,それまでに入力した全ての作業を,実行可能なPerlスクリプトとして保存してくれます.このログは作業ディレクトリのG-languageSessionLog-<日付>.plというファイルに保存されます.

   
   G > makelog
   

G-language GAEには200以上にも及ぶ豊富な解析プログラムが実装されています。さらに,全てのプログラムには複数のオプションが存在し,幅広い解析を自由にカスタマイズしながら行うことができます.しかし,当然ながらこれら200個全ての名前や使い方を覚えるのは不可能です.そこで,G-languageシェルにはマニュアルを閲覧・検索するためのhelpコマンドが実装されています。help -sのオプションでG-language GAEとBioPerlの全ドキュメントを,help -gでG-language GAEの全ドキュメントを,help -bでBioPerlの全ドキュメントをキーワード検索することができます。ここでは例としてこれから解析で使うGC skewに関連するプログラムを検索してみましょう。

   
   G >  help -g GC skew
            Found keyword "GC skew" in the following manual pages (4 hits).
   
    find_ori_ter     -   predict the replication origin and terminus in bacterial genomes
    gcsi             -   GC Skew Index: an index for strand-specific mutational bias
    gcskew           -   calculate the GC skew of the given genome
    genomicskew      -   calculate the GC skew in different regions of the given genome

マニュアルを閲覧するには,helpコマンドに直接そのプログラム名をあたえます.

   G > help genomicskew
   
    Name: genomicskew   -   calculate the GC skew in different  regions of the given genome
   
    Description:
      This program graphs the GC skew for the whole genome, coding regions,
      intergenic regions, and the third codon.
   
    Usage:
       1 = genomicskew($genome);
   
    Options:
      -divide      window number to divide into (default: 250)
      -at          1 when observing AT skew instead of GC skew (default: 0)
      -output      f for file output in directory "data",
                   g for graph output in directory "graph",
                   show for graph output and display (default: "show")
      -filename    output filename (default: "genomicskew.png" for -output=>"g",
                                             "genomicskew.csv" for -output=>"f")
      -application application to open png image (default: "gimv")
   
    Author:
       Kazuharu Arakawa (gaou@sfc.keio.ac.jp)
   
    History:
      20040610-01 updated to handle introns and exons
      20040601-01 bug fix for output=>"f" option
      20010905-01 updated options
      20010727-01 initial posting
   

他にもシェルからはPubMed(pubmedコマンド)で論文を検索したり,Entrez(entrezコマンド)からデータを検索したりすることができます。以下に例を示しますが,詳しくはhelpコマンドを使って調べてみてください.

   
   G > pubmed G-language
       3 entries found in PUBMED: (Showing up to 10 hits)
   
        1. PMID:   17572364
           Tamaki S et al. "Restauro-g: a rapid genome re-annotation system for  
           comparative genomics."
           Genomics Proteomics Bioinformatics  2007 Feb;5(1):53-8
   
        2. PMID:   16789913
           Arakawa K et al. "GPAC: benchmarking the sensitivity of genome informatics 
           analysis to genome annotation completeness."
           In Silico Biol  2006;6(1-2):49-60
   
        3. PMID:   12538262
           Arakawa K et al. "G-language Genome Analysis Environment: a workbench for 
           nucleotide sequence data mining."
           Bioinformatics  2003 Jan 22;19(2):305-6
   
   G > entrez genome Carsonella rudii
       5 entries found in GENOME: (Showing up to 10 hits)
   
        1. Accession Number:   NC_009495
           Clostridium botulinum A str. ATCC 3502, complete genome
   
        2. Accession Number:   NC_004328
           Plasmodium falciparum 3D7 chromosome 7, whole genome shotgun sequence
   
        3. Accession Number:   NC_004325
           Plasmodium falciparum 3D7 chromosome 1, complete sequence
   
        4. Accession Number:   NC_008512
           Candidatus Carsonella ruddii PV, complete genome
   
        5. Accession Number:   NC_002952
           Staphylococcus aureus subsp. aureus MRSA252, complete genome
   
   G > $gb = new G("NC_008512")
   Retrieving sequence from REFSEQ...
   
   
   Accession Number: NC_008512
   
     Length of Sequence :    159662
              A Content :     66734 (41.80%)
              T Content :     66481 (41.64%)
              G Content :     12946 (8.11%)
              C Content :     13501 (8.46%)
                 Others :         0 (0.00%)
             AT Content :    83.44%
             GC Content :    16.56%
   
   
   G > $gb->output("carsonella.gbk")
   
   G > head carsonella.gbk
   LOCUS       NC_008512  159662bp     DNA   circular  BCT      22-OCT-200
   DEFINITION  Candidatus Carsonella ruddii PV, complete genome.
   ACCESSION   NC_008512
   VERSION     NC_008512.1  GI:116334902
   KEYWORDS    .
   SOURCE      Candidatus Carsonella ruddii PV.
     ORGANISM  Candidatus Carsonella ruddii PV Bacteria; Proteobacteria;           
   Gammaproteobacteria; Candidatus
               Carsonella.
   REFERENCE   1
     AUTHORS   Nakabachi,A., Yamashita,A., Toh,H., Ishikawa,H., Dunbar,H.E.,

2.2.2 バクテリア解析システムをコマンドラインで使う

GUIで使ったバクテリア解析システムは,GCFに定義された解析を行う機能を持っていましたが,繰り返して一連の作業を行いたい場合など,GUIよりもコマンドラインが適している場合があります.そういった場合には,シェルのGコマンドにGCFファイルを引数として与えるだけで動作します.

   
   knoppix@Knoppix:~$ G ~/.glang/default.gcf
   

同様に,拡張子が".pl"のPerlスクリプトを引数に与えると,perl -MGのようにG-language GAEのモジュールをロードしたPerlコマンドと同様に扱います.

   
   knoppix@Knoppix:~$ G test.pl
   

解析プログラムを個別に使うには,関数の形でシェル内で記述します.たいていの場合は第一引数にnew G()で読み込んだゲノムファイルを与えます.

   
   G >$mgen = new G("/home/knoppix/book/chapter3/mgen.gbk");Accession Number: NC_000908
   Length of Sequence :    580074
            A Content :    200543 (34.57%)
            T Content :    195695 (33.74%)
            G Content :     92312 (15.91%)
            C Content :     91524 (15.78%)
               Others :         0 (0.00%)
           AT Content :    68.31%          
           GC Content :    31.69%
   G >gcwin($mgen);
   

オプションは,オプション名に続けて⇒で値を指定します.

   
   G > gcwin($mgen, -window=>1000);
   

3. G-languageによるバクテリアゲノム解析

   コラム: GC skewとは 
   GC skew とは,一本鎖DNA分子におけるG含量とC含量のバイアスを表す指標で,(Cの個数 - Gの個数)/(Cの
   個数 + Gの個数) の式で表されます.ワトソンとクリックによるDNA分子構造発見のきっかけともなった有名
   なシャルガフの法則は,「細胞内のAとT,そしてGとCの量はほぼ等しい」というものですが,後にこの塩基分
   配法則には第二項が加えられており,「一本鎖DNA分子においてもAとT,そしてGとCの量はほぼ等しい」と経
   験則から定められています.ただし,これはゲノム全体における現象で,局所的な領域ではAとT,そしてGとC
   の量比に偏り(バイアス)が見られます.さらに,バクテリアのうち真性細菌の多くの種ではゲノム中で明確
   にその傾向が逆転する個所がみられ,その個所が複製開始・終結地点と一致することが多いことが知られてい
   ます.GC skew という現象が発生する原因には主に二説あり,それらはリーディング鎖とラギング鎖の異なる
   突然変異確率,そしてコドン使用による変異のバイアスなどによる,と言われています.
   
   参考文献:
   
     * Lobry, J.R., 1996. Asymmetric substitution patterns in the two DNA strands of 
   bacteria. Mol. Biol. Evol. 13, 660-665 

3.1 GC skewと複製開始・終結点の関係

GC skewという現象は,バクテリアのコンプリートゲノムが初めて読み取られてから比較的早い段階で知られ,今でもゲノム全体における傾向をつかむ手法として,基礎的な解析の一つになっています.G-language GAEは,ゲノム全体を使った解析を行う上で常に考慮すべきであるGC skew関連の情報を得るツールを豊富に備えています.

それでは実際にGC skewの解析をはじめてみましょう.様々なゲノムの傾向を比較するために,今回は5つの生物のコンプリートゲノムを使用します.

生物名 種名 アクセッション番号 説明
大腸菌 Escherichia coli NC_000913 最も良く知られたモデル生物.
枯草菌 Bacillus subtilis NC_000964 納豆菌の親戚.特に強いGC skewを示す.
マイコプラズマ菌 Mycoplasma genitalium NC_000908 コンプリートゲノムが読まれている中で,最小ゲノムを持つ微生物.
シアノバクテリア Synechococcus sp. NC_005070 光合成をする細菌.
パイロコッカス菌 Pyrococcus furiosus NC_003413 超高熱耐性の古細菌.

まずはG-languageシェルを起動し,これら5つのゲノムを使用可能な状態にしましょう.これら5つの生物種のゲノムファイルはシステムにサンプルデータとして同梱されています.コンピュータのメモリ(RAM)に制限がある場合には,5つ全部を読み込むのが難しい可能性がありますので,その場合には一つずつ読み込んで解析をするようにしてください.

      
   knoppix@Knoppix:~$G
            __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/
                  G-language  Genome Analysis Environment v.1.7.2
   
                           http://www.g-language.org/
            Please cite:
                Arakawa K. et al. (2003) Bioinformatics.
                Arakawa K. et al. (2006) Journal of Pestice Science.
            License: GNU General Public License
            Copyright (C) 2001-2007 G-language Project
            Institute for Advanced Biosciences, Keio University, JAPAN
           __/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/__/           
   G >$eco = new G("ecoli", "no msg");
   G >$bsub = new G("bsub", "no msg");
   G >$mgen = new G("mgen", "no msg");
   G >$cyano = new G("cyano", "no msg");
   G >$pyro = new G("pyro", "no msg");
   G >
   

以上のコマンドを入力することで,5つのゲノムを以後の解析で利用できる形で$eco, $bsub, $mgen, $cyano, $pyroの中に格納できます.通常 new G()コマンドはゲノムデータを読み込むと同時に,全体の塩基組成の割合を表示しますが,"no msg"オプションを与えてやることでこれを表示しないように設定できます.

それでは実際にこれらのゲノムのGC skewを見てみましょう.GC skewを計算する関数はgcskew($genome)です.

   
   G >gcskew($eco)
   G >gcskew($bsub)
   G >gcskew($mgen)
   G >gcskew($cyano)
   G >gcskew($pyro)
   

以上のように入力すると,5つの新規ウインドウが立ち上がり,それぞれにグラフが表示されるはずです.GC skew関数は,ゲノムを10,000bpの領域に分割し,それぞれの領域における(C-G)/(C+G)を計算してプロットします.よって,ある領域にCがGよりも多ければ正の値をとり,逆の場合には負の値をとります.

                            大腸菌                                                         枯草菌                                                   マイコプラズマ菌

                       シアノバクテリア                                    パイロコッカス菌

大腸菌と枯草菌では相対的にGが多い領域とCが多い領域が,非常に綺麗に分かれているのが見て取れます.逆に,シアノバクテリアとパイロコッカス菌ではそういった傾向は見えません.一方,マイコプラズマ菌のグラフはデータポイント数が他のものと比べてかなり少なくなっています.どうやらマイコプラズマ菌の58万bpというゲノム長に対し,10,000bpのウインドウは大きすぎたようです(他のゲノムは2.5Mbp~4.6Mbp).オプションでウインドウサイズを1,000bpに縮小して再度グラフしてみましょう.

   
   G > gcskew($mgen, -window=>1000)
   

これで他のグラフと同じくらいの細かさで描画されました.しかし,マイコプラズマ菌のGC skewも,シアノバクテリアやパイロコッカス菌同様に明確なシフトポイントが見えません.これらの生物にはGC skewは存在しないのでしょうか?GC skewのように正と負の値をとる連続した数字の切り替わりを見るには,累積プロットが有効であることが知られています.ウインドウ毎に個別の値をプロットするのではなく,次々に値を累積して足し合わせていくのです.このように,累積GC skewを計算するには-cumulative⇒1オプションを指定します。先ほどと同様に,このオプションをつけて(マイコプラズマ菌ではウインドウ幅を1,000bpに変えて)累積GC skewを見てみましょう.

   
   G >gcskew($eco, -cumulative=>1)
   G >gcskew($bsub, -cumulative=>1)
   G >gcskew($mgen, -cumulative=>1, -window=>1000)
   G >gcskew($cyano, -cumulative=>1)
   G >gcskew($pyro, -cumulative=>1)
   

                            大腸菌                                                         枯草菌                                                   マイコプラズマ菌

                       シアノバクテリア                                    パイロコッカス菌

GC skewが明確に見えた大腸菌と枯草菌では,累積をとることによってさらにシフトポイントを強調させることができました.また,ノイズが大きくてGC skewそのものがはっきりとは観察できなかったマイコプラズマ菌やシアノバクテリアに関しても,累積GC skewを見るとシフトポイントが,そしてGC skewが弱いながらも存在していることがわかります.一方,古細菌であるパイコッカス菌は累積をとっても明確なシフトポイントは依然として現れてきません.

一般的に,真性細菌のリーディング鎖はラギング鎖に比べてGに富んでいます.よって,GC skewのシフトポイントはリーディング鎖とラギング鎖の境目,つまり,複製開始点と終結点に対応していると言われています.累積GC skewでは最高点が予測される複製開始点に,最低点が予測される複製終結点に対応しています.この二つの頂点を探すプログラムは,find_ori_ter($genome)という関数として提供されています.複製開始点と終結点が知られている大腸菌と枯草菌を使って,これらを予測してみましょう.

   
   G > find_ori_ter($eco)
   find_ori_ter:
      Window size = 1000
      Predicted Origin:   3923500
      Predicted Terminus: 1549500
   G > find_ori_ter($bsub)
   find_ori_ter:
      Window size = 1000
      Predicted Origin:   4213500
      Predicted Terminus: 1941500
   

大腸菌の複製開始点は3,923,500bp付近で,終結点は1,588,800bp付近,そして枯草菌の複製開始点は1bp付近,終結点は2,017,000bp付近です.必ずしも全てが一致しているわけではないですが,GC skewを用いた複製開始点及び終結点の予測は比較的高精度であると言えるでしょう.(枯草菌の予測された複製開始点はゲノムの最後部分を指していて,環状ゲノムなのでこれは最初の部分と等しいく,よって実際の複製開始点である1bpに極めて近いと言えます.)マイコプラズマ菌も複製開始点は1bpにあることが知られており,複製終結点は未知ですが,少なくとも累積GC skewで頂点が1bp付近にあることは現象が良く一致しています.一方で,古細菌の多くは複数の複製開始点を持つことから,パイロコッカス菌では明確なシフトポインやGC skewが見いだせなかったことも納得できます.

それでは次に,GC skewがどの部位に対する突然変異によって起きている現象なのかを調べてみましょう.GC skewが比較的強かった大腸菌と枯草菌において,ゲノム全体,コード領域のみ,遺伝子間領域のみ,そして変異が起きやすいコドンの3番目の塩基のみのGC skewを同時にプロットしてみます.この時,それぞれの領域の長さが変わってしまいますので,それぞれ全体の長さの1/250に区切った領域を一つのウインドウとします.これら全てを行ってくれるのがgenomicskew($genome)という関数です.

   
   G >genomicskew($eco)
   G >genomicskew($bsub)
   

                                                                大腸菌

                                                                枯草菌

遺伝子間領域やコドンの3番目の塩基はコード領域内及びゲノム全体と比べて長さがかなり短くなってしまいますので,その分振れ幅が大きくなっています.しかし,いずれにしてもこれらどの領域もGC skewがあることを示しています.よって,GC skewを作り出す突然変異は,ゲノム上特定部位に作用するというよりは,リーディング鎖全体にほぼ等しく起こる可能性が示唆されます.

さて,GC skewの傾向が強い種とそうでない種の違いは何でしょうか.はっきりしたことはまだ分かっていませんが,少なくとも細菌の増殖速度と何らかの関係がありそうです.枯草菌の一回の細胞分裂にかかる時間は約0.45/h,大腸菌は約0.35h,マイコプラズマ菌は約12h,そしてシアノバクテリアは約8hですが,この速度とGC skewの強度には相関がうかがえます.

   参考文献:
   
   * Rocha, E.P.C., 2004. Codon usage bias from tRNA's point of view: Redundancy, 
   specialization, and efficient decoding for translation optimization. Genome Res. 
   14, 2279-2286(増殖速度など) 
   * Freeman, J.M., et al. 1998. Patterns of Genome Organization in Bacteria. Science. 
   279, 1827a(複製開始・終結点など) 
   * Grigoriev, A., 1998. Analyzing genomes with cumulative skew diagrams. Nucleic 
   Acids Res. 26, 2286-2290(累積グラフなど) 
   * Arakawa K., et al. 2007. Noise-reduction filtering for accurate detection of 
   replication termini in bacterial genomes. FEBS Letters. 581(2), 253-258 (フーリエ変換
   を用いた最新の遺伝子開始・終結点予測法) 

3.2 シグナルオリゴ配列の傾向

   コラム: χ配列とは 
   大腸菌では一回の複製過程に数回以上,何らかの形で複製フォークの進行が妨げられ,それによりDNAの二重鎖
   切断が生じてしまいます.この時,複製フォークを修復し再度複製を続行させるために,RecBCD酵素によって
   促進される相同組換えが非常に重要な役割を演じます.RecBCDはゲノム中のχ(カイ)配列と呼ばれる8塩基 
   配列,5'-GCTGGTGG-3'を組換えのホットスポットとして認識してその働きを調節しますが,χ配列は特定の方
   向を向いている時しか作用しないため,ゲノム中で有意に複製方向を向いて存在しているシグナル配列として
   有名です.
   
   参考文献:
   
   * Kowalczykowski, S.C., 2000. Initiation of genetic recombination and recombination-
   dependent replication.Trends Biochem. Sci. 25, 156-165 
   * Uno, R. et al., 2000. The orientation bias of Chi sequence is a general tendency 
   of G-rich oligomers. Gene. 259, 207-215 

ここでは,特定のオリゴ配列が複製方向に向かっている割合,すなわち,リーディング鎖上に存在している割合を計算する方法を例示します.

リーディング鎖上に存在している割合をもとめるのですから,リーディング鎖の配列が必要です.G-language GAEでは以下のようにleading_strand($genome) 関数を使うことで複製開始点から複製終結点までの,右回りと左回りそれぞれのリーディング鎖配列を得ることができます.これらを繋げてリーディング鎖のみのゲノム配列を作ります.

   
   G > ($seq1, $seq2) = leading_strand($eco)
   G > $leading = $seq1 . $seq2
   

new G($genome)関数を使ってゲノムデータベースを読み込むと,そのゲノム塩基配列の簡単な統計が表示されますが,実はこれはseqinfo($genome)という関数を呼び出して表示しています.seqinfo($genome)を用いて,GC skewで示されるようにリーディング鎖にはG含有量が多いことを確認してみましょう.G-language GAEの関数のうち,ゲノム情報や配列を引数とする関数は自動的に型推論を行うので,ここではゲノム情報のアノテーションを含む全てを持つであるnew G($genome)から得られる構造体だけでなく,塩基配列のみを持つスカラー値を入れます.

   
   G >seqinfo($leading)
     Length of Sequence :   4639675
              A Content :   1137315 (24.51%)
              T Content :   1145883 (24.70%)
              G Content :   1216052 (26.21%)
              C Content :   1140425 (24.58%)
                 Others :         0 (0.00%)
             AT Content :    49.21%
             GC Content :    50.79%
   

ゲノム全体の傾向と比較すると,AとTの量は変化がないものの,Gは増加し,Cが低下しているのが見て取れます.χ配列は8文字中5文字がGという偏った塩基含有量を持つので,ラギング鎖よりもリーディング鎖上に多く存在していることは全体の塩基含量から推測できます.それでは,実際にχ配列の数を数えてみましょう.ある長い配列中に,短い配列が存在する個数を数えるにはoligomer_counter($genome, "oligomer")を使います.

複製方向を向いている数を数えるには,GCTGGTGGという配列がリーディング鎖中に存在している数を数えます.一方,逆方向を向いている数を数えるには,ラギング鎖中に同じ配列が存在している数を数えれば良いのですが,これはリーディング鎖中のGCTGGTGGの相補配列,CCACCAGCの数を数えることでも可能です.complement($seq)関数を使えば簡単に相補配列を得ることができます.実際に複製方向を向く割合を計算すると以下のようになります.

   
   G > print oligomer_counter($leading, "gctggtgg")761
   G > print oligomer_counter($leading, complement("gctggtgg"))247
   G > print 761+2471008G > print 761/1008 * 10075.4960317460317
   

同様のことを,任意の8塩基配列を使って計算してみましょう.いくつか試してみると,すぐに以下のことに気がつきます. 多くの8塩基配列はゲノム中に100~300個程度存在し,1008個というのは極端に多い 75%という複製方向を向く割合も同様にかなり高い 実際にこれを統計的に見るとどうなるでしょうか.まず,ゲノム中に存在する数を計算するには,ある8塩基の組み合わせが任意の配列である確率は,塩基が4通り(それぞれが25%の確率)なので0.258=0.000015程度になります.大腸菌ゲノムは4.6Mbpなので,この中に特定の8塩基オリゴ配列が一本鎖に存在している数は約70個.相補鎖に存在している数を考慮すれば計約140個.この確率から考えてもχ配列は有意に多そうだと言えます.

また,複製方向に関しては,塩基含量から配列がリーディング鎖に存在する確率とラギング鎖に存在する確率をそれぞれ

   
   F1 = [A]a[T]t[G]g[C]c
   

   
   F2= [A]a'[T]t'[G]g'[C]c'
   

([A][T][G][C]はリーディング鎖の各塩基含量,atgcはオリゴ配列中の塩基数,a't'g'c'はオリゴ配列相補配列中の塩基数) として求め,

   
   F1/(F1+F2)
   

を計算することから割合を算出できます.χ配列についてこれを計算すると,

   
   G >$f1 = 0.2451 ** 0 * 0.2470 ** 2 * 0.2621 ** 5 * 0.2458 ** 1
   G >$f2 = 0.2451 ** 2 * 0.2470 ** 0 * 0.2621 ** 1 * 0.2458 ** 5
   G >print $f1/($f1 + $f2) * 10056.7651511506961
   

と約57%の割合が予測され,75%はこれを大きく上回っていることがわかります.ただし,この解析では各塩基の含有量のみからの予測であり,実際はゲノム上の配列には複数の塩基の並びが関与するマルコフ性が知られていますので,そういったより高次の情報を使わずに結論しないよう注意が必要です.

   参考文献:
   
   * Lobry, L.R. et al., 2003. Polarisation of prokaryotic chromosomes. Curr. Opin. 
   Microbiol. 6, 101-108 
   * Arakawa K. et al., 2007. Significance of the genomic properties of Chi sites 
   validated from the distribution of all octamers in Escherichia coli. Gene. 392(1-
   2), 239-246 

3.3 全オリゴの複製方向バイアス

「シグナルオリゴ配列の傾向」ではχ配列という一つの8塩基配列について調べてみました.しかし,統計的な比較をするためには65536通りの全ての8塩基配列に対し同様の解析を行い,それら全体の傾向を見ることが望ましいでしょう.このような単純な繰り返し作業は,簡単なプログラミングを行うことで自動化することができます.G-language GAEはPerlのモジュールとしても扱えるので,実際にPerlのスクリプトを書いて全オリゴ配列の複製方向バイアスを計算してみましょう.

まず,シェルを終了する前に大腸菌のゲノムデータを保存しておきます.ゲノムデータのを保存するには,new G()で作り出した構造体からoutput()関数を呼び出し,ファイル名を引数として与えます.この時,ファイル名の拡張子(.gbk, .fatsta, .emblなど)でG-language GAEは出力形式を自動判別して適切なフォーマットで保存してくれます.ここではGenBank形式でファイルを保存しましょう.

   
   G > $eco->output("ecoli.gbk")
   

さて,それでは実際にPerlスクリプトの作成を始めましょう.Perlのスクリプト,例えば test1.pl というファイルを作成し,emacsなどのエディタで開きましょう.そして,G-language GAEを使う,ということを宣言するために,use G; と書き入れます.そして,先ほど保存した大腸菌のゲノムデータをnew G()で読み込みます.

     use G;
     $eco = new G("ecoli.gbk", "no msg");

次に,先ほどシェルで行ったように,リーディング鎖の配列を求めます.

     ($seq1, $seq2) = leading_strand($eco);
     $leading = $seq1 . $seq2;

ここからPerlによる処理が始まります.まず,リーディング鎖中に存在する全8塩基の数を数えます.これはリーディング鎖の配列の先頭から,8文字のウインドウを一文字ずつずらしながらリーディング鎖の配列の末端までどんな8塩基配列が存在しているかを数えていくことで出来ます.このような繰り返しには for 文を使います.0番目(先頭)から,8塩基のウインドウをとれるリーディング鎖の最後の部分まで,一文字ずつ,というのをfor文で表現すると以下のようになります.

     for ($i = 0; $i <= length($leading) - 8; $i ++){
     }

このfor文の中は,一文字ずつウインドウをずらす度に繰り返されます.ここでは,$i番目から8文字をリーディング鎖の配列から切り出し,その配列のカウント数を1増やす,という作業を記述すればよいですね.特定の文字列を切り出すにはsubstr()というPerlの標準関数を使い,文字列を数えるためにはハッシュを使います.

     for ($i = 0; $i <= length($leading) - 8; $i ++){
         $oligo{substr($leading, $i, 8)} ++;
     }

さて,これで全8塩基配列の,リーディング鎖中の数を数えて,%oligoというハッシュに格納することができました.あとはそれぞれの8塩基配列の総数(相補鎖の数を含む)と,複製方向に向く割合を計算すれば完成です.既にシェルでχ配列を使って行ったように,この部分は以下のように記述できます.

         $total = $oligo{$nuc} + $oligo{complement($nuc)};
         $percent = $oligo{$nuc} / $total * 100;
         printf "$nuc, $total, $percent\n";

あとはこの処理を全8塩基配列に関して行えば完成です.ハッシュの全要素に関して繰り返して行うためには,foreachを使います.

     foreach $nuc (sort keys %oligo){
         $total = $oligo{$nuc} + $oligo{complement($nuc)};
         $percent = $oligo{$nuc} / $total * 100;
         printf "$nuc, $total, $percent\n";
     }

以上でスクリプトは完成です.以下にスクリプト全体をまとめます.

     use G;
     $eco = new G("ecoli.gbk", "no msg");
     ($seq1, $seq2) = leading_strand($eco);
     $leading = $seq1 . $seq2;
 
     for ($i = 0; $i <= length($leading) - 8; $i ++){
         $oligo{substr($leading, $i, 8)} ++;
     }
 
     foreach $nuc (sort keys %oligo){
         $total = $oligo{$nuc} + $oligo{complement($nuc)};
         $percent = $oligo{$nuc} / $total * 100;
         printf "$nuc, $total, $percent\n";
     }

スクリプトを実行し,結果をout.csvというファイルに出力しましょう.カンマで区切ったデータは,CSVという形式にしておくとそのままExcelなどの表計算プログラムで開くことができて便利です.

   
   unix:~$ perl test1.pl >out.csv
   

プログラムが無事終了したら,out.csvをlessコマンドなどで見ると,以下のように結果を見ることができます.表計算ソフトなどでこのデータを読み込み,複製方向を向く割合でソートし,一番偏りが強い配列は何か,χ配列は全オリゴ配列は何番目に偏っているか,数で見るとχ配列は何番目に多いか,などを調べてみましょう.

   
   aaaaaaaa, 242, 38.4297520661157
   aaaaaaac, 351, 38.4615384615385
   aaaaaaag, 318, 47.1698113207547
   aaaaaaat, 502, 48.804780876494
   aaaaaaca, 418, 48.5645933014354
   aaaaaacc, 375, 39.4666666666667
   aaaaaacg, 387, 53.2299741602067
   aaaaaact, 286, 49.6503496503497
   aaaaaaga, 429, 51.2820512820513
   aaaaaagc, 563, 49.911190053286
   aaaaaagg, 336, 52.3809523809524
   aaaaaagt, 285, 51.9298245614035
   aaaaaata, 492, 49.7967479674797...
   

3.4 遺伝子の複製方向バイアス

ここまでに,塩基とオリゴ配列と複製方向の関係について見てきましたが,遺伝子の向きはどうなっているのでしょうか.

全オリゴ配列に関して繰り返し処理を行う時に,一塩基ずつ動かすのにfor文を,オリゴ配列毎にハッシュの中身を見るのにforeach文を使ってきましたが,ゲノム中の遺伝子を一つずつ見るために,G-language GAEでは以下の構文を使います.

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

$genome→cds()にはゲノム中の全遺伝子を指すポインタ,$cdsが格納されており,$cdsを使うことで遺伝子の様々な情報を得ることができます.例えば,主な󂂣??は以下の通りです.

表記法 説明
$genome→{$cds}→{start} 遺伝子開始位置
$genome→{$cds}→{end} 遺伝子終結位置
$genome→{$cds}→{direction} 遺伝子のある鎖(direct or complement)
$genome→{$cds}→{gene} 遺伝子名
$genome→{$cds}→{locus_tag} 遺伝子番号
$genome→{$cds}→{product} 遺伝子産物
$genome→{$cds}→{db_xref} 外部データベースのリンク
$genome→{$cds}→{translation} 翻訳されたアミノ酸配列
$genome→get_geneseq($cds) 遺伝子の塩基配列
$genome→startcodon($cds) スタートコドン
$genome→stopcodon($cds) ストップコドン
$genome→before_startcodon($cds, 100) 遺伝子上流100bpの塩基
$genome→after_startcodon($cds, 100) 遺伝子下流100bpの塩基
$genome→get_intron($cds) イントロンの配列
$genome→get_exon($cds) エキソンの配列

それでは,それぞれの遺伝子がどちらの鎖にのっているかを調べ,もしリーディング鎖にのっているのであればカウントを増やす,という方法で遺伝子が複製方向を向いている割合を計算してみましょう.遺伝子の総数はscalar($genome→cds())で求められます.どちらの鎖にあるか,を調べるためにはquery_strand($genome, <position>, -direction⇒<direction>)関数を使います.ここで,<position>は数字で表されるゲノム中の位置を指しますので,今回は遺伝子の開始位置を使います.また,<direction>には遺伝子の向きを使います.

以上を踏まえてスクリプトを作成すると,以下のようになります.

     use G;
     $eco = new G("ecoli.gbk", "no msg");
     foreach $cds ($eco->cds()){
         my $strand = query_strand($eco, $eco->{$cds}->{start},  
                                              -direction=>$eco->{$cds}->{direction});
         if($strand eq 'leading'){
             $i ++;
         }
      }
      print $i/scalar($eco->cds());

実際にこのスクリプトを実行して,遺伝子が複製方向を向いている割合を算出してみましょう.

3.5 遺伝子発現量と複製方向バイアス

   コラム: Codon Adaptation Indexとは 
   
   CAIは,リボソーム関連遺伝子などの高い発現遺伝子の同義コドン使用頻度の偏りをもとにして,これら高発現
   遺伝子のコドンバイアスと観察対象遺伝子のコドンバイアスの相関から相対的な発現量を推定する手法です.
   当然CAIのみで定量的かつ正確な遺伝子発現量を予測することは不可能ですが,ある程度発現量の目安を得るこ
   とができ,少なくともコドンバイアスが,方向性を持って強く存在している真性細菌においては,その目安と
   しての有用性は一般的に認められています.
   
   参考文献:
   
   * Sharp, P.M. et al., 1987. The codon Adaptation Index--a measure of directional 
   synonymous codon usage bias, and its potential applications. Nucleic Acids Res. 15, 
   1281-1295 
   

さて,これまでにゲノム中の様々な構成部分と複製方向の関係を見てきましたが,最後に,より機能的な情報である遺伝子の発現量に注目して,複製方向と遺伝子発現量の関係を解析してみましょう.遺伝子の発現量を見るにはマイクロアレーなどのハイスループットな網羅的実験手法によるデータを活用するのが理想的ですが,ここではゲノム情報から遺伝子発現量を予測する指標として代表的なCodon Adaptation Index (CAI)を用います.

ゲノム全体の傾向は,少なからず個々の遺伝子や短い塩基配列に影響を与えます.それはGC skewなどの塩基組成の偏りが存在していたり,リーディング鎖とラギング鎖における遺伝子の数が異なったりするからです.大腸菌はリッチな環境では約30分に一回細胞分裂しますし,複製という機構はバクテリアゲノムにおいて非常に大きな意味を持ちます.ここでは,リーディング鎖とラギング鎖によって遺伝子の発現量にも何かしらの影響がある,という仮説を検証してみます.

G-language GAEでCAIを使う場合,cai($genome)関数を使います.CAIは各遺伝子毎に0から1の間の値(0が高発現遺伝子と無相関で,1が相関)で得られます.cai($genome)を使うと全遺伝子のCAI値のリストが得られるほか,$genome→{$cds}→{cai}に新たに遺伝子のCAI値が追加されます.今回はこの値を使用し,リストは不要なため,-outputオプションでリストを出力しないように調整します.

     cai($genome, -output=>"/dev/null");

さて,ゲノム領域中の傾向を見るには,ゲノム全体を細かい領域に分けて,それぞれにおける値をグラフとしてプロットするのが有効であることはこれまでのGC skewなどの解析でみてきた通りです.また,傾向をより顕著にするためには累積プロットが適しています.以上を踏まえて,CAI値をプロットしてみましょう.しかし,ここで二点留意する必要があります.

まず,CAI値は0から1の値をとりますが,累積プロットを行う為には正と負の値を持つものでなければならないので,これを補正してやる必要があります.大腸菌の全遺伝子のCAI値の平均は約0.485なので,実際のCAI値からこの値を引いてやることで,平均が0,最大値が約0.5,最小値が約-0.5となります.

次に,今まではGC含量などだけだったのでDNAの一本鎖だけ考慮すればよかったのですが,遺伝子はどちらの鎖にも存在しますので,プロットをリーディング鎖とラギング鎖の二つに分ける必要があります.どちらの鎖を使うか判別する部分には,既に遺伝子の方向性の解析で使ったquery_strand($genome)関数が使えます.

各遺伝子で,平均値0に補正したCAI値を計算し,リーディング鎖とラギング鎖に分けます.これを後にグラフにするので,それぞれをアレーを使って格納します.アレーに要素を追加するにはPerl標準関数のpush()を使います.CAI値だけでなく,遺伝子のゲノム中の位置も保存しておきましょう.以上をプログラムにすると以下のようになります.

     foreach $cds ($eco->cds()){
         my $strand = query_strand($eco, $eco->{$cds}->{start}, 
                                              -direction=>$eco->{$cds}->{direction});
         if ($strand eq 'leading'){
             $cai_leading += $eco->{$cds}->{cai} - 0.485;
         }else{
                   $cai_lagging += $eco->{$cds}->{cai} - 0.485;
         }
         push(@data1, $cai_leading);
         push(@data2, $cai_lagging);
         push(@pos, $eco->{$cds}->{start});
     }

最後に,得られたデータをグラフにします.G-language GAEでは簡単にグラフを作成することができるgrapher()という関数があるので,これに値を入れて行きます.それぞれ配列の参照の形で,各データのX座標のリスト(上で格納したゲノム中の位置)と,各データのY座標のリスト(上で格納したリーディング鎖のCAI値とラギング鎖のCAI値)を入れるだけです.

     grapher(\@pos, \@data1, \@data2);

これだけで作業ディレクトリ内部のgraphというフォルダ内にgraph.pngという名前でグラフの画像が作成されます.グラフをさらにわかりやすくするように,タイトル,データラベル,X軸やY軸のラベルを追加します.

     grapher(\@pos, \@data1, \@data2, -x1=>"leading", -x2=>"lagging",
                                         -title=>"caiskew", -x=>"position", -y=>"cai");

以上をまとめると以下のようなプログラムになります.

     use G;
     $eco = new G("ecoli.gbk", "no msg");
     cai($eco, -output=>"/dev/null");
 
     foreach $cds ($eco->cds()){
         my $strand = query_strand($eco, $eco->{$cds}->{start}, 
                                               -direction=>$eco->{$cds}->{direction});
         if ($strand eq 'leading'){
             $cai_leading += $eco->{$cds}->{cai} - 0.485;
         }else{
                   $cai_lagging += $eco->{$cds}->{cai} - 0.485;
         }
         push(@data1, $cai_leading);
         push(@data2, $cai_lagging);
         push(@pos, $eco->{$cds}->{start});
     }
 
      grapher(\@pos, \@data1, \@data2, -x1=>"leading", -x2=>"lagging",
                                        -title=>"caiskew", -x=>"position", -y=>"cai");

グラフを表示するには,displayコマンドを使います.

   
   unix:~$ display graph/graph.png
   

すると,以下のようなグラフが表示されます.

既に見てきたように,大腸菌の複製開始点は3,923,500bp付近で,終結点は1,588,800bp付近です.まず,ラギング鎖を見ると,ほぼ全域が負の値になっています.これはラギング鎖全体を通じて平均以下のCAI値を持つ遺伝子が多いこと,つまり,ラギング鎖中の遺伝子の発現量が相対的に少ないことを示しています.

一方,リーディング鎖のグラフは大きく変動しています.そして,複製開始点と終結点との関連を見ると,複製終結点周辺ではグラフが下降を続け,複製開始点付近では大幅に上昇していっているのがわかります.よって,リーディング鎖では複製開始点付近に高発現遺伝子が集中しており,徐々に終結点に向かうにつれて発現量は低くなっていっているようです.

さて,このグラフでは,3.5Mbp付近でリーディング鎖の値が急激に上昇しています.CAI値はリボソーム関連遺伝子を高発現遺伝子と定義していて,リボソーム関連遺伝子はかたまって遺伝子クラスターを形成していることが多いので,その影響がでていそうです.そこで,リボソーム関連遺伝子を除いて同じグラフを作ってみましょう.foreach文の一番最初に,

     next if ($eco->{$cds}->{product} =~ /ribosom/);

と記入するだけで,リボソーム関連の遺伝子は除去できます.グラフはどのようにかわりましたか?

   参考文献:
   
   * Daubin, V. et al., 2003. G+C3 Structuring Along the Genome: A Common Feature in
   Prokaryotes. Molecular Biology and Evolution. 20, 471-483
   
   * Arakawa K. et al., 2007. Selection effects on the positioning of genes and gene
   structures from the interplay of replication and transcription in bacterial 
   genomes. Evolutionary Bioinformatics. 3, 279-286

4. おわりに

本章では,G-language GAEとPerlを使って,実際にプログラミングなどもしながらバクテリアゲノムの解析を行ってきました.プログラミングに慣れていなければ難しいところもあったかもしれません.しかし,逆に言えば以上に示した解析例は,とくに後半はかなり複雑な解析もありましたが,いずれも20行未満のスクリプトです.実際に研究を行う場合はもっとプログラミングが必要になるでしょうが,いずれにしてもG-language GAEではこのように複雑な解析を効率よく,比較的少ない手続きで行うことを可能にします.

それでは,あとはこのG-language Genome Analysis Environmentを使ってあなたが生命の真理を探究してみてください.

introductiontog-languagejapanese.txt · Last modified: 2014/01/18 07:44 (external edit)