(57)【要約】      （修正有）【課題】DNAアレイ等の遺伝子発現プロファイリング実験結果から得られた類似の発現パターンを示す遺伝子の上流領域に共通して存在する、遺伝子発現制御配列のモチーフを推定する方法を提供する。【解決手段】（１）遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子の上流領域を入力配列とし、入力配列を１０〜４０ｂｐの長さの構成要素配列に分解し、（２）各構成要素配列に関し、構成要素配列に相同な配列を検索し、検索された相同な配列及び構成要素配列からなる配列群に関してモチーフを作成し、（３）作成されたモチーフから配列の保存性を表わす統計量を算出し、（４）ゲノム配列からサンプリングした、モチーフと同一の数及び長さの配列からなる配列群から（３）の工程と同じ方法で算出された統計量と、（３）の工程で算出された統計量との有意差を検定し、有意差のあるものを、遺伝子発現制御配列の推定モチーフとして選択する。【選択図】図１
2007108949.tif 000002
【特許請求の範囲】
【請求項１】
下記の工程を有することを特徴とする、遺伝子発現プロファイリングの結果から遺伝子発現制御配列のモチーフを推定する方法。
（１）遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子の上流領域を入力配列とし、入力配列を１０〜４０ｂｐの長さの構成要素配列に分解する。
（２）各構成要素配列に関し下記の工程を行う。
（２−１）入力配列における、構成要素配列に相同な配列を検索する。
（２−２）検索された相同な配列及び構成要素配列からなる配列群に関してモチーフを作成する。
（３）作成されたモチーフから配列の保存性を表わす統計量を算出する。
（４）ゲノム配列からランダムにサンプリングした、モチーフと同一の数及び長さの配列からなる配列群から（３）の工程と同じ方法で算出された統計量と、（３）の工程で算出された統計量との有意差を検定し、有意差のあるものを、推定される遺伝子発現制御配列のモチーフとして選択する。
【請求項２】
（２−１）の工程で検索された相同な配列の数に基づき構成要素配列を選択する請求項１記載の方法。
【請求項３】
統計量が、塩基配列間の一致度、相同性又はアラインメント長である請求項１又は２記載の方法。
【請求項４】
統計量が、塩基配列を構成する塩基の情報量又は出現確率である請求項１又は２記載の方法。
【請求項５】
遺伝子発現プロファイリングに用いる遺伝子が原核生物のものである請求項１〜４のいずれか１項に記載の方法。
【請求項６】
原核生物が、アミノ酸生産能及び／又は有機酸生産能を保持する微生物である請求項５記載の方法。
【請求項７】
遺伝子発現プロファイリングが、DNAアレイ、EST解析、及び、RT−PCRから選ばれるものである、請求項１〜６のいずれか１項に記載の方法。
【請求項８】
遺伝子発現制御配列が、プロモーター配列、転写制御因子結合配列、small RNA、及び、トランスポゾンのいずれかである、請求項１〜７のいずれか１項に記載の方法。
【請求項９】
請求項１〜８のいずれか１項の方法で遺伝子発現制御配列のモチーフを推定し、推定された遺伝子発現制御配列のモチーフを含む複数のモチーフについて距離行列を作成し、作成された距離行列に基づいてモチーフを分類することを含む、モチーフを分類する方法であって、モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出し、算出された相同性又は距離に基づいてモチーフの距離行列を作成することを特徴とする前記方法。
【請求項１０】
請求項１〜８のいずれか１項の方法で遺伝子発現制御配列のモチーフを推定し、推定された遺伝子発現制御配列のモチーフと機能既知のモチーフとの相同性又は距離を求め、相同性のある又は距離の近い機能既知のモチーフの機能を、機能を推定するモチーフの機能と予測することを含む、モチーフの機能を予測する方法であって、モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出することにより相同性を求めることを特徴とする前記方法。
【請求項１１】
請求項１〜８のいずれかの方法により遺伝子発現制御配列のモチーフを推定し、推定された遺伝子発現制御配列のモチーフまたは請求項１０に記載の方法により機能が予測されたモチーフに類似する配列の検索方法であって、
モチーフを学習データとして隠れマルコフモデルを作成し、
作成された隠れマルコフモデルを用い、遺伝子の上流領域及びゲノム配列について検索を行い、
ゲノム配列から検索された配列の有するスコアに対して有意に大きいスコアを有する上流領域から検索された配列を、類似する配列として選択する
ことを含む前記方法。
【請求項１２】
請求項１０に記載の方法によりモチーフの機能を予測し、機能が予測されたモチーフを上流領域に持つ遺伝子からなる遺伝子群を、共通の遺伝子発現制御機構を持つものとして結び付けることを含む、転写ネットワークの作成方法。
【請求項１３】
遺伝子群が、前記モチーフと類似した配列を上流領域に持ち、かつ遺伝子発現プロファイリングで発現が変化したことが確認された遺伝子をさらに含み、前記モチーフと類似した配列は請求項１１に記載の方法で検索されたものである、請求項１２に記載の方法。
【請求項１４】
下記の処理をコンピューターに実行させることを特徴とする、遺伝子発現プロファイリングの結果から遺伝子発現制御配列のモチーフを推定するためのプログラム。
（１）遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子の上流領域を入力配列とし、入力配列を１０〜４０ｂｐの長さの構成要素配列に分解する。
（２）各構成要素配列に関し下記の工程を行う。
（２−１）入力配列における、構成要素配列に相同な配列を検索する。
（２−２）検索された相同な配列及び構成要素配列からなる配列群に関してモチーフを作成する。
（３）作成されたモチーフから配列の保存性を表わす統計量を算出する。
（４）ゲノム配列からランダムにサンプリングした、モチーフと同一の数及び長さの配列からなる配列群から（３）の工程と同じ方法で算出された統計量と、（３）の工程で算出された統計量との有意差を検定し、有意差のあるものを、推定される遺伝子発現制御配列のモチーフとして選択する。
【請求項１５】
モチーフ間の相同性又は距離を求める方法であって、モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出することを含む前記方法。
【請求項１６】
複数のモチーフについて距離行列を作成し、作成された距離行列に基づいてモチーフを分類することを含む、モチーフを分類する方法であって、請求項１５の方法でモチーフ間の相同性又は距離を求め、求めた相同性又は距離に基づいてモチーフの距離行列を作成することを特徴とする前記方法。
【請求項１７】
機能を推定するモチーフと機能既知のモチーフとの相同性又は距離を求め、相同性のある又は距離の近い機能既知のモチーフの機能を、機能を推定するモチーフの機能と予測することを含む、モチーフの機能を予測する方法であって、請求項１５の方法でモチーフ間の相同性又は距離を求めることを特徴とする前記方法。
【請求項１８】
モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出する処理をコンピューターに実行させることを特徴とする、モチーフ間の相同性又は距離を求めるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、DNAアレイ等の遺伝子発現プロファイリング実験結果から得られた類似の発現パターンを示す遺伝子の上流領域に共通して存在する、遺伝子発現制御配列を推定する方法、ならびに、それを利用して、モチーフの分類、モチーフの機能予測、モチーフ類似配列の検索、及び、転写ネットワーク構造の推定をする方法、ならびに、モチーフ間の相同性又は距離を求める方法、ならびに、それらの方法をコンピューターに実行させるためのプログラムに関する。
【背景技術】
【０００２】
生体内で遺伝子の発現は種々の因子の影響を受けることが知られている。そして、生物にとっての影響の大きさなどに依存して、遺伝子発現の関係は階層構造をとる。近年、DNA アレイ等の技術の発展により、種々の環境における遺伝子発現プロファイルを調べることで、遺伝子発現関係を推定することができるようになってきている。
【０００３】
遺伝子は転写制御因子により発現制御を受けることが知られている。転写制御因子は特定のDNA配列に対し強い親和性を持って結合する。転写制御因子により認識される複数のDNA配列から、その転写制御因子に特徴的なパターン配列あるいはコンセンサス配列を検出することが出来る。転写制御因子、及び、特徴的なパターン配列を同定することは、産業上有用であることが知られている遺伝子の発現制御に直接関わることから、遺伝子工学上有用な知見となる。
【０００４】
また、DNAアレイを用いた実験により、特定の条件下で類似の発現パターンを示す複数の遺伝子を検出することが可能である。類似の発現パターンを示す遺伝子群は、共通した転写制御因子により支配されていることが推定される。従って、そのような遺伝子群の発現を支配する転写制御因子を同定することが要求されるが、現在のところ、DNAアレイ実験の結果から直接、転写制御因子を同定する技術は存在しない。そこで、この様な遺伝子の上流領域から転写制御因子結合配列を発見し、その配列を基にバイオインフォマティクス及び、実験的な手法を活用し、遺伝子発現制御に関わる転写制御因子を同定する手段が提案されている（非特許文献１〜３、特許文献１〜２）。
【０００５】
転写制御因子が結合する配列は、通常10-30mer程度の長さで、高度に保存されている。このことから、類似の発現パターンを示す遺伝子群の上流配列から、局所的に保存されたパターンを発見することで、転写制御因子が結合する配列を推定することが出来ると考えられている。この様なパターン発見の目的に利用されるアルゴリズムとして、最尤法やギブスサンプリングといったものが代表的であり、これらのアルゴリズムを実装した代表的なソフトウエアとしてMEME (Multiple EM for Motif Elicitation)(非特許文献５), Gibbs Sampler(非特許文献６)などが存在する。
【特許文献１】WO2004/053106
【特許文献２】WO01/86577
【非特許文献１】Mol. Microbiol. 48, 699-712, 2003
【非特許文献２】Genome Res. 10: 744-757
【非特許文献３】J. Mol. Biol. 284, 241-254
【非特許文献４】HYPERLINK "http://www.genome.org/content/vol12/issue10/" Genome Res. 12: 1523-1532
【非特許文献５】Proceedings of the Second International Conference on Intelligent Systems for Molecular Biology, pp. 28-36, AAAI Press, Menlo Park, California, 1994
【非特許文献６】Science, 262, 208-214
【発明の開示】
【発明が解決しようとする課題】
【０００６】
既存技術を用いた場合、遺伝子の発現関係を推定するためには数多くのDNAアレイ実験結果が要求される。その結果、試験研究に莫大なコストがかかる上に、推定された発現関係は恒常的に観察される関係に限定される。そのため、ある特定の環境でのみ出現する、条件特異的な遺伝子発現関係を推定するのは困難である。
【０００７】
また、与えられた配列中に存在する全てのパターンを、既存技術を単独で実行することにより検出することはアルゴリズム上不向きである。また、与えられた配列の長さが増加するにつれて、パターン検出の感度、精度が悪化する。例えば、Gibbs Samplerソフトウエアを利用する際は発見するべきパターンの長さを指定することが要求される(Lawrence,
C.E., Altschul, S.F., Bogouski, M.S., Liu, J.S., Neuwald, A.F., and Wooten, J.C. (1993), "Detecting Subtle Sequence Signals: A Gibbs Sampling Strategy for Multiple Alignment," Science, 262, 208-214., http://melina.hgc.jp/expl_par.html#gibbs_top)。しかし、同定すべき転写制御因子が未知の場合、結合するDNA配列の長さも未知である事から、長さを指定する要求を満たすことは困難である。
【０００８】
また、発見されたDNAモチーフは、その後、DNA-蛋白質相互作用を利用した、転写制御因子を同定するための実験に利用する、あるいは、既に知られているDNAモチーフと比較し、遺伝子発現に関与する転写制御因子を推定することを前提としている。MEMEやGibbs Samplerにより発見されたDNA モチーフが生物学的に意味を持ちうるものか、現在のところ判断する方法が少ない(Natalia Poluliakh and Kenta Nakai, Extraction of Biological Motifs by Gibbs Sampler from the Promoters of Homo Sapiens, Saccharomyces cerevisiae and Bacillus subtilis Genome Informatics 14: 406-407 (2003))。発見したモチーフの検定は全体の実験精度という点で非常に重要である。
【０００９】
MEMEやGibbs Samplerに代表される従来の方法を単独で利用した場合、入力配列に対して可能性のあるすべてのDNAモチーフを検出することが出来なかった。また、これらの方法により発見されたDNAモチーフについて、入力配列における統計的有意性は示された。しかし、通常ゲノム配列と比べ、入力配列が非常に短いことから、ゲノム配列レベルでの有意性について判断する指標が存在しなかった。更に、機能既知のDNAモチーフと機能未知のDNAモチーフを比較するために、BLASTやFASTAに代表される配列比較法は、任意の比率でA,G,T,Cが存在する塩基配列集合同士を比較することが出来ないため、DNAモチーフを代表する塩基を強制的に決める必要があり、DNAモチーフ間の類似度を間接的にしか表現できなかった。これらの諸課題のため、遺伝子発現プロファイリングから遺伝子発現制御配列を推定すること、さらにはそれらを基に転写ネットワーク構造を推定することが困難であった。
【課題を解決するための手段】
【００１０】
本発明者らは、上記課題を解決するために、ある生物由来の複数のDNA配列から、長さを先に指定すること無しに、局所的に保存されたDNAモチーフを全て検出し、その生物種のゲノム配列と比較して、統計的に有意に保存されていることを検定することにより、遺伝子発現制御配列のモチーフを推定する方法を見出した。更に、DNAモチーフを、シャノンエントロピーを用いて数値行列に変換し、両端のギャップペナルティを含まない準グローバルアラインメント法により、DNAモチーフ間の類似度をベクトルの内積の最大値として直接的に計測する方法を見出した。これらの手法を活用し、発現変化の見られた遺伝子の上流に存在し、遺伝子制御配列として機能するDNAモチーフを発見するとともに、その情報を基に転写ネットワーク構造を推定する方法を見出した。
【００１１】
本発明は、下記のとおりである。
［１］  下記の工程を有することを特徴とする、遺伝子発現プロファイリングの結果から遺伝子発現制御配列のモチーフを推定する方法。
（１）遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子の上流領域を入力配列とし、入力配列を１０〜４０ｂｐの長さの構成要素配列に分解する。
（２）各構成要素配列に関し下記の工程を行う。
（２−１）入力配列における、構成要素配列に相同な配列を検索する。
（２−２）検索された相同な配列及び構成要素配列からなる配列群に関してモチーフを作成する。
（３）作成されたモチーフから配列の保存性を表わす統計量を算出する。
（４）ゲノム配列からランダムにサンプリングした、モチーフと同一の数及び長さの配列からなる配列群から（３）の工程と同じ方法で算出された統計量と、（３）の工程で算出された統計量との有意差を検定し、有意差のあるものを、推定される遺伝子発現制御配列のモチーフとして選択する。
［２］  （２−１）の工程で検索された相同な配列の数に基づき構成要素配列を選択する［１］記載の方法。
［３］  統計量が、塩基配列間の一致度、相同性又はアラインメント長である［１］又は［２］記載の方法。
［４］  統計量が、塩基配列を構成する塩基の情報量又は出現確率である［１］又は［２］記載の方法。
［５］  遺伝子発現プロファイリングに用いる遺伝子が原核生物のものである［１］〜［４］のいずれか１項に記載の方法。
［６］  原核生物が、アミノ酸生産能及び／又は有機酸生産能を保持する微生物である［５］記載の方法。
［７］  遺伝子発現プロファイリングが、DNAアレイ、EST解析、及び、RT−PCRから選ばれるものである、［１］〜［６］のいずれか１項に記載の方法。
［８］  遺伝子発現制御配列が、プロモーター配列、転写制御因子結合配列、small RNA、及び、トランスポゾンのいずれかである、［１］〜［７］のいずれか１項に記載の方法。
［９］  ［１］〜［８］のいずれか１項の方法で遺伝子発現制御配列のモチーフを推定し、推定された遺伝子発現制御配列のモチーフを含む複数のモチーフについて距離行列を作成し、作成された距離行列に基づいてモチーフを分類することを含む、モチーフを分類する方法であって、モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出し、算出された相同性又は距離に基づいてモチーフの距離行列を作成することを特徴とする前記方法。
［１０］  ［１］〜［８］のいずれか１項の方法で遺伝子発現制御配列のモチーフを推定し、推定された遺伝子発現制御配列のモチーフと機能既知のモチーフとの相同性又は距離を求め、相同性のある又は距離の近い機能既知のモチーフの機能を、機能を推定するモチーフの機能と予測することを含む、モチーフの機能を予測する方法であって、モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出することにより相同性を求めることを特徴とする前記方法。
［１１］  ［１］〜［８］のいずれかの方法により遺伝子発現制御配列のモチーフを推定し、推定された遺伝子発現制御配列のモチーフまたは［１０］に記載の方法により機能が予測されたモチーフに類似する配列の検索方法であって、
モチーフを学習データとして隠れマルコフモデルを作成し、
作成された隠れマルコフモデルを用い、遺伝子の上流領域及びゲノム配列について検索を行い、
ゲノム配列から検索された配列の有するスコアに対して有意に大きいスコアを有する上流
領域から検索された配列を、類似する配列として選択する
ことを含む前記方法。
［１２］  ［１０］に記載の方法によりモチーフの機能を予測し、機能が予測されたモチーフを上流領域に持つ遺伝子からなる遺伝子群を、共通の遺伝子発現制御機構を持つものとして結び付けることを含む、転写ネットワークの作成方法。
［１３］  遺伝子群が、前記モチーフと類似した配列を上流領域に持ち、かつ遺伝子発現プロファイリングで発現が変化したことが確認された遺伝子をさらに含み、前記モチーフと類似した配列は［１１］に記載の方法で検索されたものである、［１２］に記載の方法。
【００１２】
［１４］  下記の処理をコンピューターに実行させることを特徴とする、遺伝子発現プロファイリングの結果から遺伝子発現制御配列のモチーフを推定するためのプログラム。
（１）遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子の上流領域を入力配列とし、入力配列を１０〜４０ｂｐの長さの構成要素配列に分解する。
（２）各構成要素配列に関し下記の工程を行う。
（２−１）入力配列における、構成要素配列に相同な配列を検索する。
（２−２）検索された相同な配列及び構成要素配列からなる配列群に関してモチーフを作成する。
（３）作成されたモチーフから配列の保存性を表わす統計量を算出する。
（４）ゲノム配列からランダムにサンプリングした、モチーフと同一の数及び長さの配列からなる配列群から（３）の工程と同じ方法で算出された統計量と、（３）の工程で算出された統計量との有意差を検定し、有意差のあるものを、推定される遺伝子発現制御配列のモチーフとして選択する。
【００１３】
［１５］  モチーフ間の相同性又は距離を求める方法であって、モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出することを含む前記方法。
［１６］  複数のモチーフについて距離行列を作成し、作成された距離行列に基づいてモチーフを分類することを含む、モチーフを分類する方法であって、［１５］の方法でモチーフ間の相同性又は距離を求め、求めた相同性又は距離に基づいてモチーフの距離行列を作成することを特徴とする前記方法。
［１７］  機能を推定するモチーフと機能既知のモチーフとの相同性又は距離を求め、相同性のある又は距離の近い機能既知のモチーフの機能を、機能を推定するモチーフの機能と予測することを含む、モチーフの機能を予測する方法であって、［１５］の方法でモチーフ間の相同性又は距離を求めることを特徴とする前記方法。
【００１４】
［１８］  モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出する処理をコンピューターに実行させることを特徴とする、モチーフ間の相同性又は距離を求めるためのプログラム。
【発明の効果】
【００１５】
本方法によって、遺伝子プロファイリング解析の質を大幅に高め、現在活発に進行中のポストゲノム研究の迅速な遂行に寄与することが可能である。また、得られたDNAモチーフ情報を、DNAプローブ配列の設計、新規転写制御因子の同定と機能解明に活用することが可能であり、遺伝子工学における技術の高精度化に大きく寄与することが可能である。
【発明を実施するための最良の形態】
【００１６】
表１を用いて、以下、本発明の一実施の形態を説明する。
【００１７】
【表１】
2007108949.tif 000003
【００１８】
本発明は、通常には、処理1：構成要素分解、処理２：構成要素毎のモチーフ作成（処理２−１：局所的相同配列検索、処理２−２：配列アラインメント）、処理３：モチーフの統計量計算、処理４：モチーフに関する検定、処理５：モチーフの相同性又は距離の算出(処理５−１：モチーフの数値化、処理５−２：相同性又は距離の計算)、処理６：モチーフの分類、処理７：モチーフの機能予測、処理８：モチーフ類似配列の検索、処理９：ネットワークの生成から選ばれる処理の可能の組合せであって、少なくとも処理１〜４又は処理５を含む組合せから構成される。
【００１９】
処理１及び処理２により、入力配列から可能性のある全てのモチーフ候補を作成し、処理３により全モチーフ候補に対して統計値を求め、それを基に処理４により有意性を検定することで、モチーフを推定（抽出）する。処理５によりモチーフの相同性又は距離を算出する。処理６により、抽出されたモチーフを含む複数のモチーフに関し、モチーフを分
類する、すなわち、モチーフのコンセンサス配列パターンが類似しているものをひとまとめにする。処理７により、機能既知のモチーフとの比較により、抽出されたモチーフの機能を推定する。処理８により、抽出されたDNAモチーフ様配列を遺伝子の上流領域に持つ遺伝子を網羅的に検索する。処理９により、処理７（及び処理８）の結果を基に遺伝子ネットワークを推定する。
【００２０】
既に機能既知のDNAモチーフが必ず遺伝子ネットワークに含まれることがわかっている場合は処理６と同時に処理７を実施することが可能である。また、処理７において、既知の機能が推定されない場合もある。
【００２１】
本発明において、遺伝子発現制御配列とは、転写結合因子結合配列、プロモーター配列、small RNA、トランスポゾン、の全長およびその部分配列を指す。
【００２２】
処理１では、遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子の上流領域を入力配列とし、入力配列を１０〜４０ｂｐの長さの構成要素配列に分解する（図１）。
【００２３】
遺伝子発現プロファイリングに用いる遺伝子は、特に限定されないが、原核生物のものであることが好ましい。原核生物は、例えば、アミノ酸生産能及び／又は有機酸生産能を保持する微生物である。アミノ酸生産能及び／又は有機酸生産能を保持するとは、アミノ酸及び／又は有機酸の生産に使用できることを意味する。
【００２４】
遺伝子発現プロファイルとは、DNAアレイ、EST解析、RT-PCRなどの手法により得られる、複数の遺伝子の発現状態を示すものである。遺伝子プロファイルは、公開されたデータベースから得ることもできる。
【００２５】
公開されたデータベースとしては、奈良先端科学技術大学院大学より公開されているDNAアレイデータのデータベース（GenoBase（ecoli.aist-nara.ac.jp/gb4/search/xp_analysis/2_components/all.html））が知られている。このデータベースのDNAアレイデータはCy3, Cy5を利用した二色蛍光型のDNAアレイデータであり、コントロールとして、野生型のE. coliから抽出したRNAを用いている。
【００２６】
遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子は、遺伝子プロファイリングの結果（遺伝子プロファイル）から、発現パターンを解析する通常の方法により選択することができる。例えば、遺伝子発現プロファイルデータを、階層型クラスター化法、K-meansクラスター化法、自己組織化マップといったクラスタリング法を用いて解析し、同一、または隣り合うクラスターに分類されたものを選択する（Steen Kundsen著、塩島聡、松本治、辻本豪三監訳、わかる！使える！DNAマイクロアレイデータ解析入門、羊土社）。
【００２７】
遺伝子発現プロファイルデータ、例えばDNAアレイデータの解析には、既存のソフトウェア、例えばGene Cluster 3.0とMicrosoft Excel（登録商標）を利用することができる（de Hoon MJ, Imoto S, Nolan J, Miyano S. 2004. Open source clustering software.
Bioinformatics. 20:1453-1454.）。
【００２８】
クラスタリングの解析の前に、検定に基づく抽出を行ってもよい。例えば、野生株対野性株、野生株対目的遺伝子欠損株のDNAアレイデータを比較し、目的遺伝子欠損株で発現量が増加している遺伝子を抽出するために以下のような処理を行ってもよい。各スポット比を対数変換した後、遺伝子毎のスポット比の対数値の平均値と標準偏差を算出する。たとえば複数種類のDNAアレイデータにおいて、1種類でも有意な発現変動が見られた遺伝子
を抽出することとし、多重性を考慮した検定を行うことができる。具体的には、分散分析の手法にのっとり、遺伝子毎の不偏分散値（下記参照）を算出し、それを利用してt検定を実施するための統計値を算出する（坪井  達夫著。Excelで学ぶ統計・統計で学ぶExcel、エーアイ出版、p88,123）。更に、実験回数の少ないDNAアレイデータの場合には、t検定の結果得られたp値を一律に３〜５倍（例えば4倍）にするという補正を行い、その補正値を基に有意差検定の結果を評価してもよい。有意差を判定する際の有意水準は比較パターン数で補正してもよい(Bonferroni補正)。
【００２９】
遺伝子毎の不偏分散Ve=偏差平方和Se / 自由度Ne
偏差平方和Se=Σアレイデータ｛測定値（対数変換したスポット値）−標本平均（対数変換したスポット値の平均）｝2
自由度Ne=全アレイにおけるスポット数−実験条件数
統計量＝｜Δ標本平均｜/ √[不偏分散Ve×(1/スポット数) ×２]
有意水準＝0.05/3または0.01/3
【００３０】
更に、二色蛍光型のDNAアレイデータを解析する場合には、ほとんどの遺伝子の発現には差は無いことが前提であることを利用し、野生株対野性株のDNAアレイデータでは、二つのスポット間で差が大きすぎるものは除外し、そのほかのデータでは二つのスポット比の対数値が０とは有意に異なるもののみを採用してもよい。具体的には各DNAアレイのスポット比の対数値の平均と標準偏差を計算し、スポット比の対数値が、スポット平均±標準偏差の2倍の範囲に収まらない遺伝子を抽出し、野生株対野性株のDNAアレイデータではこれらの値を示す遺伝子は解析対象から除外する。また、そのほかのデータでは、上記検定で有意差があると判定された遺伝子群が、抽出した遺伝子群に含まれることを確認する。
【００３１】
これらの二つの検定により、目的遺伝子欠損株において有意に発現が変動した遺伝子を抽出した後、それらの遺伝子を対象として、クラスタリング解析を行うことができる。本発明では、入力配列として遺伝子発現プロファイリングにおいて類似の発現パターンを示す遺伝子の上流領域からなる複数のDNA配列を用いる。
【００３２】
類似の発現パターンを示す遺伝子から、目的遺伝子破壊株において、野生株に比べて有意に発現が変化した遺伝子の内、オペロンの先頭にある遺伝子又は、オペロンを形成していない遺伝子を特に選択してもよい。
【００３３】
上流領域は、推定しようとする発現制御配列に応じて選択すればよいが、通常には、開始コドンの上流100〜1000bpから下流の10〜200bpを選択する。例えば、開始コドン上流500bpから下流の100bpを入力配列とすることができる。
【００３４】
構成要素配列の長さは、１０〜４０ｂｐから適宜選択される。同一の入力配列を、複数の長さの構成要素配列に分解してもよい。例えば、入力配列を20mer、25mer及び30merの構成要素配列に分解してもよい。分解は、既存のソフトウェア、例えばEMBOSS wordcount
(EMBOSS: The European Molecular Biology Open Software Suite (2000) Rice,P. Longden,I. and Bleasby,A. Trends in Genetics 16, (6) pp276-277)を用いて行うことができる。
【００３５】
処理２−１では、上記のようにして得られた構成要素配列のそれぞれについて、入力配列における、構成要素配列に相同な配列を検索する（図１）。検索は、分解された構成要素配列と相同な配列を検索する他は、通常の相同性配列検索と同様でよい。この際には、入力配列に対し、構成要素及び構成要素と相同な配列が、どの位置にどの位の数が存在するかを明らかにすることが好ましい。そのために、先ず、構成要素配列と相同な配列を検
出する。相同である事の基準は、構成要素配列の長さに応じて選択すればよいが、例えば、構成要素配列（問い合わせ配列）が30merの場合異なる塩基の数が13以下、25merの場合10以下、20merの場合7以下とすることができる。本処理は、既存のソフトウェア、例えばEMBOSS fuzznucを利用して行うことができる。fuzznucの出力結果として、各問い合わせ配列について、検索対象配列中に指定した条件に合う配列の合計数とその位置情報、および配列情報が出力される。遺伝子の発現変化の大多数を支配すると推定される遺伝子制御配列を発見することが目的である場合には、検索対象配列中に条件に合う配列が一定数以下（例えば遺伝子数が２１のとき合計4本以下）の問い合わせ配列はこれ以降の処理の対象外としてもよい。
【００３６】
処理２−２では、構成要素配列毎に、検索された相同な配列及び構成要素配列からなる配列群に関してモチーフを作成する（図１）。この処理は、処理２−１で得られた構成要素配列及びその相同な配列からなる配列群を対照とする他は通常のモチーフ作成と同様でよい。通常には、上記の様にして得られた配列群についてマルチプルアラインメントを実施し、DNAモチーフを作成する。本処理は、MEME(Timothy L. Bailey and Charles Elkan,
Fitting a mixture model by expectation maximization to discover motifs in biopolymers, Proceedings of the Second International Conference on Intelligent Systems for Molecular Biology, pp. 28-36, AAAI Press, Menlo Park, California, 1994.)、clustalw、gibbs sampler等のプログラムを利用して行うことが可能である。MEMEを実施する際には、入力配列群に対し保存された配列が必ず一つ存在する、あるいは、ゼロ個もしくは一つ存在する、あるいは、複数個存在する、とする3種類の計算方法をすべて実施してもよい。モチーフは、通常、モチーフ部分の配列の集合として得られる。また、MEMEを実行する際に利用する遺伝子の塩基組成に関するマルコフモデルとして、E. coliの非遺伝子領域における６次のマルコフモデルを作成し、利用した。ここでは、非遺伝子領域はE. coliゲノムにおける蛋白質遺伝子以外の部分を指す。蛋白質遺伝子のゲノムにおける領域情報は既に公開されているゲノム情報から得た。非遺伝子領域配列を1mer-5mer長さの構成要素が、それぞれの長さにおいてどのように分布するかを調べることで非遺伝子領域の６次のマルコフモデルを作成することが出来る。
【００３７】
マルチプルアラインメントの実施により、一定の長さの問い合わせ配列群に対して、より短い配列からなる出力結果が得られることがある。このことから、DNAモチーフを構成する複数の配列中に複数の同一の配列が存在する可能性が生じる。この様な重複を確認するために、モチーフを構成する配列が、もとの入力配列のどこに位置するものであるかを確認することが好ましい。本プロセスを実施するためにEMBOSS fuzznucを利用することができる。同じモチーフに、全く同じ配列が存在し、尚且つ、同じ入力配列の同じ位置に由来するものである場合、それらは配列が重複しているものとみなすことが出来るので、一つにまとめる。
【００３８】
処理３では、作成されたモチーフから配列の保存性を表わす統計量を算出する（図２）。配列の保存性を表わす統計量の例としては、シャノンエントロピーによるInformation Content(IC)、ローカルアラインメント又はグローバルアラインメントによる配列間の一致度等が挙げられる。
【００３９】
ICとは、シャノンエントロピーを用いてDNA配列を数値に変換した結果を指す（Schneider, T.D., Stormo, G.D., Gold, L., Ehrenfeucht, A. 1986. Information content of binding sites on nucleotide sequences. J. Mol. Biol. 188: 415-431.）。
【００４０】
ローカルアラインメント又はグローバルアラインメントによる配列間の一致度等として、さらに具体的には、モチーフ候補を構成する全塩基配列間の一致度、相同性やアラインメント長等が挙げられる。ここでの一致度とは、二つの塩基配列を比較したときに、同一
の文字の出現頻度である。また、相同性とは、Smith-Watermanアルゴリズム、あるいはNeedleman-Wunschアルゴリズムにより二つの塩基配列を比較したときにPAMやBLOSUMといったスコア行列を用いて計算される、遺伝的な変異の起こりやすさを考慮した類似性である(Smith, T.F., Waterman, M.S. (1981) Identification of common molecular subsequences J. Mol. Biol. 147: 195-197、Needleman, S. B. and Wunsch, C. D. 1970. A general method applicable to the search for similarities in the amino acid sequence of two proteins. J. Mol. Biol. 48: 443-453)。また、アラインメント長とは、Smith-Watermanアルゴリズム、あるいはNeedleman-Wunschアルゴリズムにより二つの塩基配列を比較したときの、アラインメントの開始点から終了点までの長さである。
【００４１】
配列の保存性を表わす統計量の算出は、処理２で得られたモチーフを対象とする他は通常の統計量計算と同様でよい。
【００４２】
例えば、得られたモチーフはN行M列の文字行列とみなすことがきる。各列における文字の保存性はシャノンエントロピーを利用して数値化することができる。例えば、シャノンエントロピーによるInformation Content (IC)は、下記の式により算出できる（Schneider TD, Stormo GD, Gold L, Ehrenfeucht A., Information content of binding sites on
nucleotide sequences. J Mol Biol. 1986 Apr 5;188(3):415-31.）。
【００４３】
【数１】
2007108949.tif 000004
【００４４】
E. coliのゲノムにおけるA,G,T,Cの出現頻度はそれぞれ25％であるが、微生物の種類によってはそれらの出現頻度は均等ではない。その場合は、出現頻度に応じてシャノンエントロピーにおける各塩基の存在確率を変更することが可能である。シャノンエントロピーを用いて表された各残基のInformation Content (IC)の合計値を情報量と呼ぶ。
【００４５】
また、行毎の配列の保存性は、Smith-Waterman法に代表されるようなローカルアラインメント法、またはNeedleman-Walsh法に代表されるようなグローバルアラインメント法により配列間の一致度、相同性、あるいは、アラインメントされた長さを利用して表現することができる（Needleman, S. B. and Wunsch, C. D. (1970) J. Mol. Biol. 48, 443-453; Smith TF, Waterman MS (1981) J. Mol. Biol 147(1);195-7）。
【００４６】
また、得られたモチーフの出現確率に基づいてモチーフの選択をしてもよい。例えば、配列の保存性を示す統計量が同じ数値であった場合、偶然に得られる確率が低いものが何らかの生物学的意味を持つ可能性が高いと考えられるので、偶然に得られる確率が低いものを選択する。この選択は、処理４の前に行っても後に行ってもよい。
得られたモチーフの出現確率は、例えば、二項分布を用いて表現することができる。
【００４７】
【数２】
2007108949.tif 000005
【００４８】
ただし、出現確率をｐ、nをモチーフの構成配列数、mをモチーフにおけるAあるいはGあるいはTあるいはCの数、また、各塩基の事前確率を0.25とする。得られた確率値の対数を取り、マイナス１をかけることでモチーフ出現スコアとして表現することができる。
【００４９】
【数３】
2007108949.tif 000006
【００５０】
処理４では、ゲノム配列から得たランダムにサンプリングした、モチーフと同一の数及び長さの配列からなる配列群から処理３と同じ方法で算出された統計量と、処理３で算出された統計量との有意差を検定し、有意差のあるものを、推定される遺伝子発現制御配列のモチーフとして選択する。検定自体は通常の方法と同様でよい。例えば、処理３により得られたモチーフの統計値と、ゲノム配列からランダムに配列をサンプリングして作成した擬似モチーフの統計値に関して、有意差の有無をｔ検定に代表される統計検定法により判定する。具体的な検定方法としては以下の二通りの方法が挙げられる。まず第一に、被検定モチーフと同じ数の配列、同じ長さの擬似モチーフをゲノム配列から作成し、列ごとの文字の保存性、Smith-Waterman法により算出した配列間の一致度を利用した行毎の配列の保存性についてノンパラメトリック検定（Wilcoxon符号つき順位和検定）によりモチーフの統計値の有意差検定を実施する。ただ一度の擬似モチーフとの比較ではゲノム配列における文字列の偏りに起因する検定結果の誤りが起こりうるので、同様の操作を最低２００回以上繰り返し、その四分の三以上有意差があると判定されたものに関して、ゲノム配列中に有意に保存されたDNAモチーフであると判定する。更に、統計検定の結果有意と判定されたモチーフのうち1残基当たりの平均情報量が1以上であるモチーフを選択してもよい。別法として、ゲノム配列からランダムに10-40merの配列を200本以上サンプリングし、長さごとに、それら配列間の保存性を調べ、その平均値をもって、ゲノム配列由来の真の値とする（図３）。この真の値を利用することで配列間の保存性について、一変数の仮説検定を実施することが可能となり、ノンパラメトリック検定（Wilcoxon符号つき順位和検定）を採用する方法に比べ計算時間の節約が期待できる。
【００５１】
ゲノム配列は、遺伝子発現プロファイリングの対象となった遺伝子が由来する生物のゲノム配列である。
【００５２】
処理５では、モチーフをそれぞれシャノンエントロピーにより数値化し、数値化されたモチーフ間の相同性又は距離を最適アラインメント法により算出することによりモチーフ間の相同性又は距離を求める。
【００５３】
シャノンエントロピーによる数値化は、それ自体公知の方法で行うことができる（図４）。これにより、DNAモチーフ間の距離の計算に際して、DNAモチーフにおける各部位の文字の保存性が考慮されることになる。数値化されたモチーフは、通常には、行列により表わされる。
【００５４】
最適アラインメントは、シャロンエントロピーにより数値化された文字列を対象とする他は、通常の最適アラインメント法と同様に行って相同性又は距離を求めればよいが、両端のギャップペナルティを含まない準グローバルアラインメント法(図５、６)、又は、ダイナミックプログラミング法を用いることが好ましい。
【００５５】
処理６では、処理１〜４により得られたモチーフを含む複数のモチーフについて、処理５により算出された相同性又は距離に基づいて距離行列を作成し、作成された距離行列に基づいてモチーフを分類する。分類は、処理５で得られた相同性又は距離に基づいて距離行列を作成する他は、通常の距離行列に基づく分類と同様でよい。例えば、処理1から４により得られたモチーフを含む複数のモチーフ間の距離を、処理５により計算し、距離行列を作成し、それをもとにクラスタリング法による分類を実施する。DNAモチーフ間の距
離の計算に際して、DNAモチーフにおける各部位の文字の保存性を考慮するために、シャノンエントロピーを利用してDNAモチーフを表現する(図４)。DNAモチーフ間の相同性を計算するために両端のギャップペナルティを含まない準グローバルアラインメント法により計算する(図５、６)。各DNAモチーフ間の相同性に基づくスコアを算出し、その逆数を計算することでモチーフ間の距離行列を得る。また別法としてダイナミックプログラミング法によるDNAモチーフ間の距離の測定も可能であり、それを基にして距離行列を作成してもよい。処理６に供するDNAモチーフの中に、すでに機能が知られているDNAモチーフを供することも可能である。距離行列に基づくクラスタリング法は種々知られており、どのような方法を用いても、DNAモチーフは分類可能である。例えばフィンランド大学で開発されたMATLAB（商標）用SOM Toolboxに実装されている自己組織化マップ法が知られている（Neural Network Tool for Data Mining: SOM Toolbox: a paper on the computational
load and applicability of SOM Toolbox 2.0 in TOOLMET 2000:）。後記実施例のSOMによるクラスタリングの結果を図７に示す。  クラスタリング解析の結果、同じクラスターに分類された複数のDNAモチーフが出現した場合は情報量が大きく、かつ、DNAモチーフを構成する配列群が偶然に出現する確率の低いDNAモチーフを選択することにより、機能未知のモチーフの絞込みを行うことができる。
【００５６】
処理７では、機能を推定するモチーフと機能既知のモチーフとの相同性又は距離を処理５により求め、相同性のある又は距離の近い機能既知のモチーフの機能を、機能を推定するモチーフの機能と予測する。機能予測は、処理５により得られる相同性又は距離を用いる他は、通常の機能予測と同様でよい。通常には、次のとおりである。データベースとして、機能既知の複数のDNAモチーフを用意し、問い合わせとして機能未知のDNAモチーフを用意する。データベース中のモチーフと問い合わせのモチーフはともにシャノンエントロピーを利用して数値に置き換える。問い合わせモチーフとデータベースモチーフ間の相同性を両端のギャップペナルティを含まない準グローバルアラインメント法により計算し、相同性の高いモチーフを検出する。モチーフ間のアラインメントを確認後、データベース中に類似のモチーフが確認できる場合、問い合わせモチーフの機能は、データベース中の機能既知のモチーフと同様であることが推定される。
【００５７】
処理８では、処理４で選択された及び／又はさらに処理７で機能が予測されたモチーフを学習データとして隠れマルコフモデルを作成し、作成された隠れマルコフモデルを用い、遺伝子の上流領域及びゲノム配列について検索を行い、ゲノム配列から検索された配列の有するスコアに対して有意に大きいスコアを有する上流領域から検索された配列を、類似する配列として選択する。通常には次のとおりである。モチーフを構成する配列群から隠れマルコフモデル(HMM)を作成する。隠れマルコフモデルの作成はMATLAB（商標） Bioinformatics Toolbox, hmmprofstruct関数、hmmprofestimate関数を利用することにより可能である。次に、作成されたHMMを利用し、検索対象となるゲノム配列上のHMMの相同領域を、MATLAB（商標） Bioinformatics Toolbox, hmmprofalign関数を利用して、ゲノム配列を100bpごとに分割し、分割された配列すべてから検出する。この解析によりHMMに対する、与えられたゲノムにおけるアラインメントスコアの平均値と標準偏差を算出することが可能である。次に、遺伝子の上流領域上のHMMの相同領域を探索するためのウィンドウ解析(Window size=30bp, Step size=10bp)を実施する。検出された相同領域のスコアがゲノム配列から得られたスコアに対して有意に大きいスコアであることがz検定により示された場合、ゲノム配列上、有意に保存されたDNAモチーフと相同な領域であると推定することができる（図９）。
【００５８】
処理９では、処理７で機能が予測されたモチーフ、及び、任意に、処理８で検索されたモチーフ類似配列を上流領域に持つ遺伝子からなる遺伝子群を、共通の遺伝子発現制御機構を持つものと結び付けて、転写ネットワークを作成する。転写ネットワークの作成は、処理７及び任意に処理８による結果を用いる他は通常のネットワークの作成と同様でよい
。通常には次のとおりである。処理７によりDNAアレイデータの発現プロファイルと関係するDNAモチーフを上流領域に持つ遺伝子群、及び、処理８でDNAモチーフと極めてよく似た配列を上流領域に持ち、DNAアレイデータにより発現が変化したことが確認された遺伝子群が共通の遺伝子発現制御機構を持つものとして線で結ばれる。
【００５９】
後記実施例で、推測された転写制御配列及びそれに類似した配列（図９）を例にして、遺伝子ネットワークを予測すると図１０に示すような結果となる。
【００６０】
また、本発明は、上記の処理をコンピューターに実行させるためのプログラムを提供する。上記の処理をコンピューターに実行させるためのプログラムは、通常のプログラム化の方法に従って作成することができる。
【００６１】
本発明のプログラムがコンピューターに実行させる処理には、キーボード等の入力装置や、他の装置からのデータ受信装置から入力する処理、ディスプレイなどの出力装置や、他の装置へのデータ送信装置へ出力する処理などの入出力処理など、プログラム化における通常の付加的な処理が含まれていてもよい。
【００６２】
また、本発明にかかるプログラムを、コンピューター読み取り可能な記録媒体に格納することもできる。ここで、この「記録媒体」は、フロッピー（登録商標）ディスク、光磁気ディスク、ＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等の任意の「可搬用の物理媒体」や、各種コンピュータシステムに内蔵されるＲＯＭ、ＲＡＭ、ＨＤ等の任意の「固定用の物理媒体」、あるいは、ＬＡＮ、ＷＡＮ、インターネットに代表されるネットワークを介してプログラムを送信する場合の通信回線や搬送波のように、短期にプログラムを保持する「通信媒体」を含むものとする。
【００６３】
また、「プログラム」は、任意の言語や記述方法にて記述されたデータ処理方法であり、ソースコードやバイナリコード等の形式を問わない。なお、「プログラム」は必ずしも単一的に構成されるものに限られず、複数のモジュールやライブラリとして分散構成されるものや、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）に代表される別個のプログラムと協働してその機能を達成するものをも含む。なお、実施の形態に示した各装置において記録媒体を読みとるための具体的な構成、読み取り手順、あるいは、読み取り後のインストール手順等については、周知の構成や手順を用いることができる。
【００６４】
以下、本発明を実施例によりさらに説明する。
【実施例１】
【００６５】
大腸菌の野生株、arcA遺伝子破壊株、及びarcB遺伝子破壊株におけるDNAアレイデータにおいて、野生株に比べ遺伝子破壊株において有意に発現が増加した遺伝子の上流配列からのarcA結合配列モチーフ発見、及び、それらを利用した遺伝子ネットワーク予測を行った。
【００６６】
DNAアレイデータは奈良先端大学より公開されているGenoBaseより入手した（http://ecoli.aist-nara.ac.jp/gb4/search/xp_analysis/2_components/all.html）。本DNA アレイデータはCy3, Cy5を利用した二色蛍光型のDNAアレイデータであり、コントロールとして、野生型のE. coliから抽出したRNAを用いている。DNAアレイデータの解析にはGene Cluster 3.0とMicrosoft Excel（登録商標）を利用した（de Hoon MJ, Imoto S, Nolan J, Miyano S. 2004. Open source clustering software. Bioinformatics. 20:1453-1454.）。野生株対野性株、野生株対ArcA欠損株、野生株対ArcB欠損株のDNAアレイデータを比較し、ArcA欠損株及びArcB欠損株でともに発現量が増加している遺伝子を抽出するために以下のような処理を行った。各スポット比を対数変換した後、遺伝子毎のスポット比の対数値
の平均値と標準偏差を算出した。3種類のDNAアレイデータにおいて、1種類でも有意な発現変動が見られた遺伝子を抽出することとし、多重性を考慮した検定を行うこととした。具体的には、分散分析の手法にのっとり、遺伝子毎の不偏分散値（下記参照）を算出し、それを利用してt検定を実施するための統計値を算出した（坪井  達夫著。Excelで学ぶ統計・統計で学ぶExcel、エーアイ出版、p88,123）。更に、使用したDNAアレイデータでは実験回数が少ないことを考慮し、t検定の結果得られたp値を一律に4倍にするという補正を行い、その補正値を基に有意差検定の結果を評価した。有意差を判定する際の有意水準は比較パターン数で補正した(Bonferroni補正)。
【００６７】
主な計算式は下記のとおりである。
遺伝子毎の不偏分散Ve=偏差平方和Se / 自由度Ne
偏差平方和Se=Σアレイデータ｛測定値（対数変換したスポット値）−標本平均（対数変換したスポット値の平均）｝2
自由度Ne=全アレイにおけるスポット数−実験条件数
統計量＝｜Δ標本平均｜/ √[不偏分散Ve×(1/スポット数) ×２]
有意水準＝0.05/3または0.01/3
【００６８】
更に、二色蛍光型のDNAアレイデータを解析する際、ほとんどの遺伝子の発現には差は無いことが前提である。このことを利用し、野生株対野性株のDNAアレイデータでは、二つのスポット間で差が大きすぎるものは除外し、そのほかのデータでは二つのスポット比の対数値が０とは有意に異なるもののみを採用することとした。具体的には各DNAアレイのスポット比の対数値の平均と標準偏差を計算した。そして、スポット比の対数値が、スポット平均±標準偏差の2倍の範囲に収まらない遺伝子を抽出した。野生株対野性株のDNAアレイデータではこれらの値を示す遺伝子は解析対象から除外した。また、そのほかのデータでは、上記検定で有意差があると判定された遺伝子群が、抽出した遺伝子群に含まれることを確認した。
【００６９】
これらの二つの検定により、ArcA欠損株、あるいはArcB欠損株、もしくはそれら両方において有意に発現が変動した遺伝子を抽出した後、それらの遺伝子を対象として、クラスタリング解析を行った。本実施例では、K-meansクラスター解析法により遺伝子の分類を行い、ArcA欠損株、ArcB欠損株でともに遺伝子発現量が有意に増加した遺伝子のみを選択した。
【００７０】
従来の方法でMEMEは最も成功したソフトウエアとして知られている。そこで、arcA遺伝子破壊株、及び、arcB遺伝子破壊株においてともに発現が上昇した遺伝子で、オペロンの先頭にある遺伝子、または、オペロンを形成していない遺伝子（acnB, acs, bglJ, fepA,
glcD, gltA, hdeA, icdA, lpdA, osmB, purE, sdhC, sucA, ydcI, ydgT, yejG, yhhX, yojH）の開始コドン上流500bpから下流側100bpを入力配列として、MEME ver. 3によるDNAモチーフの推定を実行した。結果を図１１に示す。図１１中、(A)は、全ての入力配列に１個のArcA結合モチーフが存在すると仮定した場合、(B)は、各入力配列に０個または１個のArcA結合モチーフが存在すると仮定した場合、(C)は、各入力配列に０個以上のArcA結合モチーフが存在すると仮定した場合の結果である。また、(D)は、McGuireらにより提案されたArcA結合モチーフ（McGuire, A.M., Hughes, J.D., Church, G.M. 2000. Conservation of DNA regulatory motifs and discovery of new motifs in microbial genomes, Genome Res. 10: 744-757.）である。
【００７１】
いずれの条件においても、これまで知られているArcA結合モチーフと類似したDNAモチーフは発見できなかった。
【００７２】
同じ入力配列を用いて、本発明の方法に従って、DNAモチーフの推定を行った。
まず、入力配列を20mer、25mer及び30merの構成要素に分解した（処理１−１）。本プロセスを実施するためにEMBOSS wordcountを用いた。この結果、20merの構成要素が10398種類、25merの構成要素が10308種類、30merの構成要素が10218種類得られた。
【００７３】
続いて入力配列に対し、構成要素及び構成要素と相同な配列が、どの位置にどの位の数が存在するかを明らかにした(処理１−２)。そのために、先ず、構成要素と相同な配列を検出した。相同である事の基準は、異なる塩基の数が、構成要素（問い合わせ配列）が30merの場合13以下、25merの場合10以下、20merの場合7以下とした。本プロセスを実施するためにEMBOSS fuzznucを利用した。fuzznucの出力結果として、各問い合わせ配列について、検索対象配列中に指定した条件に合う配列の合計数とその位置情報、および配列情報が出力される。本実施例では、21遺伝子の発現変化の大多数を支配する遺伝子制御配列を発見することが目的であるので、検索対象配列中に条件に合う配列が合計4本以下の問い合わせ配列はこれ以降の処理の対象外とした。この結果、20merの構成要素10398個のうち317個の構成要素が、25merの構成要素10308個のうち374個の構成要素が、30merの構成要素10218個のうち375個の構成要素が検索対象配列中に、条件に合う配列を合計5本以上持つことがわかった。
【００７４】
上記の様にして得られた相同な配列群についてマルチプルアラインメントを実施し、DNAモチーフを作成した（処理２）。本プロセスを実施するためにMEME(Timothy L. Bailey and Charles Elkan, Fitting a mixture model by expectation maximization to discover motifs in biopolymers, Proceedings of the Second International Conference on Intelligent Systems for Molecular Biology, pp. 28-36, AAAI Press, Menlo Park, California, 1994.)を用いた。MEMEを実施する際には、入力配列群に対し保存された配列が必ず一つ存在する、あるいは、ゼロ個もしくは一つ存在する、あるいは、複数個存在する、とする3種類の計算方法をすべて実施した。この結果、処理1で得られた合計1066個の構成要素に対してその3倍の数である3198個のMEME出力結果が得られた。
【００７５】
MEMEの実施により20mer,25mer,30merの問い合わせ配列群に対して、より短い配列からなる出力結果が得られることがある。このことから、DNAモチーフを構成する複数の配列中に複数の同一の配列が存在する可能性が生じる。この様な重複を確認するために、モチーフを構成する配列が、もとの入力配列のどこに位置するものであるかを確認した。本プロセスを実施するためにEMBOSS fuzznuc(EMBOSS: The European Molecular Biology Open
Software Suite (2000) Rice,P. Longden,I. and Bleasby,A. Trends in Genetics 16, (6) pp276-277)を利用した。同じモチーフに、全く同じ配列が存在し、尚且つ、同じ入力配列の同じ位置に由来するものである場合、それらは配列が重複しているものとみなすことが出来るので、一つにまとめた。
【００７６】
得られたモチーフはN行M列の文字行列とみなすことができる。そこで、得られたモチーフについて、統計量として、列毎の文字の保存性、及び、行毎の配列の保存性を算出した（処理３）。
【００７７】
各列における文字の保存性については、シャノンエントロピーによるInformation Content (IC)を、下記の式により算出した（Schneider TD, Stormo GD, Gold L, Ehrenfeucht
A., Information content of binding sites on nucleotide sequences. J Mol Biol. 1986 Apr 5;188(3):415-31.）。
【００７８】
【数４】
2007108949.tif 000007
【００７９】
E. coliのゲノムにおけるA,G,T,Cの出現頻度はそれぞれ25％とした。シャノンエントロピーを用いて表された各残基のICの合計値をモチーフの統計量とした。
【００８０】
また、行毎の配列の保存性は、Smith-Waterman法により配列間の一致度を算出した（Needleman, S. B. and Wunsch, C. D. (1970) J. Mol. Biol. 48, 443-453; Smith TF, Waterman MS (1981) J. Mol. Biol 147(1);195-7）。
【００８１】
得られたモチーフの統計量と、ゲノム配列からランダムに配列をサンプリングして作成した擬似モチーフの統計量に関して、有意差の有無を統計検定法により判定した（処理４）。具体的には、被検定モチーフと同じ数の配列、同じ長さの擬似モチーフをゲノム配列から作成し、列ごとの文字の保存性、及び、行毎の配列の保存性についてノンパラメトリック検定（Wilcoxon符号つき順位和検定）によりモチーフの統計値の有意差検定を実施した。ただ一度の擬似モチーフとの比較ではゲノム配列における文字列の偏りに起因する検定結果の誤りが起こりうるので、同様の操作を最低２００回以上繰り返し、１５０回以上有意差があると判定されたものに関して、ゲノム配列中に有意に保存されたDNAモチーフであると判定した。本検定により合計129個のモチーフが有意なモチーフと判定された。
【００８２】
以上の結果から、従来の方法であるMEMEによるDNAモチーフの推定では検索できないDNAモチーフであっても、本発明の方法によれば検索できることが分かる。
【実施例２】
【００８３】
実施例１で統計検定の結果有意と判定されたモチーフのうち1残基当たりの平均情報量が1以上である20個のモチーフを更なる処理の対象とした。
【００８４】
（１）モチーフの分類
処理対象の複数のDNAモチーフ間の距離を計算し、距離行列を作成し、それをもとにクラスタリング法による分類を実施した（処理６）。DNAモチーフ間の距離の計算に際して、DNAモチーフにおける各部位の文字の保存性を考慮するために、シャノンエントロピーを利用してDNAモチーフを表現した（図４）。DNAモチーフ間の相同性を計算するためには両端のギャップペナルティを含まない準グローバルアラインメント法により計算した(図５、６)。各DNAモチーフ間の相同性に基づくスコアを算出し、その逆数を計算することでモチーフ間の距離行列を得た。距離行列に基づくクラスタリング法として、フィンランド大学で開発されたMATLAB（商標）用SOM Toolboxに実装されている自己組織化マップ法（Neural Network Tool for Data Mining: SOM Toolbox: a paper on the computational load and applicability of SOM Toolbox 2.0 in TOOLMET 2000）により、DNAモチーフを分類した。本実施例では上記の２０個のモチーフに加え、機能が既知であるArcA結合モチーフ、Mlc結合モチーフを加えてSOMによるクラスタリングを実施した。結果を図７に示す。
【００８５】
（２）モチーフ類似配列の検索
また、処理対象の複数のDNAモチーフを学習データとして、隠れマルコフモデル(HMM)をMATLAB（商標） Bioinformatics Toolboxの hmmprofstruct関数、hmmprofestimate関数を利用して作成した（処理７）。次に、作成されたHMMを利用し、検索対象となるゲノム配列（すなわち、E. coliゲノム配列）上のHMMの相同領域を、MATLAB（商標） Bioinformatics Toolboxのhmmprofalign関数を利用して、ゲノム配列を100bpごとに分割し、分割され
た配列すべてから検出した。この解析によりHMMに対する、与えられたゲノムにおけるアラインメントスコアの平均値と標準偏差を算出した。次に、遺伝子の上流領域上のHMMの相同領域を探索するためのウィンドウ解析(Window size=30bp, Step size=10bp)を実施した。検出された相同領域のスコアがゲノム配列から得られたスコアに対して有意に大きいスコアであることがz検定により示された場合、ゲノム配列上、有意に保存されたDNAモチーフと相同な領域であると推定した。この結果、glcB遺伝子上流にArcA結合モチーフと有意に類似した配列を発見した。
【実施例３】
【００８６】
実施例１により得られたDNAモチーフの生物学的な有意性を示すため、機知の知見との比較を実施した。ArcA結合モチーフを推定する実施例１では、ArcA破壊株で発現が有意に上昇し、且つ、ArcAにより発現が制御されていることが知られている遺伝子（acnB, gltA, lpdA,sdhCB, sucCD, glcB, icd）の内、glcBとicd遺伝子上流からはArcA結合モチーフが検出できなかったが、それ以外の既知の知見を反映したモチーフが推測された(図９)。
【００８７】
以上の点を考慮し、実施例２の（１）で分類された、DNAアレイデータの発現プロファイルと関係するDNAモチーフを上流領域に持つ遺伝子群、及び、実施例２の（２）で検索された、DNAモチーフと極めてよく似た配列を上流領域に持ち、DNAアレイデータにより発現が変化したことが確認された遺伝子群が共通の遺伝子発現制御機構を持つものとして線で結ぶことにより、遺伝子ネットワークを予測すると図１０に示すような結果となる。
【００８８】
この結果、既知の知見を反映すると同時に、新規のArcA支配下の遺伝子(ydcI, yejG, yojH, b1163)を推定することが出来た。この予測を利用し、新規の育種対象となる遺伝子候補を抽出することが可能となる。また、ArcA欠損は発酵生産に有利な形質と不利な形質を同時に与えることが示唆されている。本手法によりArcA結合配列が予測されるので、プロモーター改変技術を応用し、ArcA結合配列を改変することで、特定の遺伝子のみをArcAの制御から外し、ArcA欠損の発酵生産に有利な形質のみを選択的に付与することが可能になる。
【図面の簡単な説明】
【００８９】
【図１】処理１〜２の実行手順の説明図。
【図２】処理３の実行手順の説明図。
【図３】DNAモチーフの統計検定方法の説明図。
【図４】文字集合のシャノンエントロピーによる数値化の説明図。
【図５】曖昧な文字列を比較する方法（その１）の説明図。
【図６】曖昧な文字列を比較する方法（その２）の説明図。
【図７】SOMによるDNAモチーフの分類結果。
【図８】HMMのスコア分布。
【図９】DNAモチーフの既知の知見との合致の検討結果。
【図１０】ArcAによる制御モデルのネットワーク図。
【図１１】MEMEを単独で利用した場合のDNAモチーフ抽出結果と既知ArcA結合モチーフ。(A)全ての入力配列に１個のArcA結合モチーフが存在すると仮定した場合の結果。(B)各入力配列に０個または１個のArcA結合モチーフが存在すると仮定した場合の結果。(C)各入力配列に０個以上のArcA結合モチーフが存在すると仮定した場合の結果。(D)既知ArcA結合モチーフ。
【図１】
2007108949.tif 000008
【図２】
2007108949.tif 000009
【図３】
2007108949.tif 000010
【図４】
2007108949.tif 000011
【図５】
2007108949.tif 000012
【図６】
2007108949.tif 000013
【図７】
2007108949.tif 000014
【図８】
2007108949.tif 000015
【図９】
2007108949.tif 000016
【図１０】
2007108949.tif 000017
【図１１】
2007108949.tif 000018
【配列表】
2007108949000001.app