スコア行列を用いた転写因子結合部位の予測

概要

転写因子結合部位の予測と探索は, 転写調節機構の解明の上で主要な研究対象の1つである. 当課題では, これらの問題をコンピュータを用いて解析する手法として【スコア行列を用いたプロモーター領域からの転写因子結合部位の探索】を行う。これを解くプログラムを実装し、出芽酵母のゲノム配列データを用いて、その有効性を確認する。

入力として、転写因子8種類に対する結合部位モチーフ配列群と、8つの遺伝子のプロモーター配列を与える。これより頻度行列およびスコア行列を作成し、転写因子結合部位を予測せよ。なお、バックグラウンドとして、実際のS. cerevisiaeゲノムから算出した以下の塩基出現頻度を用いる。

塩基	頻度
A	7519429
C	4637676
G	4637676
T	7519429

プログラムの作成ができたら、閾値を統計学的に検討するとよい。（発展課題）

実験ガイド

プログラミング環境

本講義では、開発環境としてreplit使用する。プログラム本体の提出もreplitを利用する。ここを参照してreplit内のKeioUniv/3jikkenteamに参加。各課題のプログラムを作成したのちに、結果の出力をまとめたファイルを(‘result.txt’などとして)同じプロジェクト上に配置し、右上のsubmitボタンで提出。提出後も期限内であれば変更可能である。

配列データとサンプルプログラムのダウンロード

ここから、配列データとサンプルプログラムを全部まとめてダウンロードし、展開する。中身の説明は下の通り。

プロジェクトkadai2にmotifフォルダ、seqフォルダをドラッグ&ドロップしてインポートする。（sample-prog_kai.cを利用してプログラムを作る場合はこれも同じプロジェクトにインポートする。）

【配列データ】

motifフォルダに結合部位配列のファイル8個が、seqフォルダにプロモーター配列ファイルが入っています。

説明	ファイル名
8つの遺伝子のプロモーター配列	`seq/promoters`
転写因子 MATa1 の結合部位配列	`motif/MATa1`
転写因子 MATalpha2 の結合部位配列	`motif/MATalpha2`
転写因子 MCM1 の結合部位配列	`motif/MCM1`
転写因子 MIG1 の結合部位配列	`motif/MIG1`
転写因子 PHO4 の結合部位配列	`motif/PHO4`
転写因子 RCS1 の結合部位配列	`motif/RCS1`
転写因子 ROX1 の結合部位配列	`motif/ROX1`
転写因子 TAF の結合部位配列	`motif/TAF`

【サンプルプログラム(3ファイル)】

sample-prog_kai.c
本課題では、プログラムを動かすごとに異なるファイル（転写因子結合部位配列ファイル）を指定することになるため、毎回プログラムを書き換えるのではなく、引数からファイル名を得る仕様にすると良い。
sample-prog_kai.cは結合部位配列ファイルおよびプロモーター配列ファイルを引数で指定して、内容を格納するプログラムです。引数として打ち込んだ2つの文字列が変数として渡され、1つ目の内容を motif[]に、2つ目の内容を pro[] に格納します。使いたい人は使ってください。
引数指定および実行の仕方は「プログラム上のヒント」の「読み込むファイルを引数で指定する」を参照してください。以下の2つの.cファイルも同様。
readargs-print.c
引数に入力された文字列（*argv[]に格納）を表示するプログラムです。argcやargvの扱い方がわかりやすいと思います。
readargs-store.c
引数に入力されたファイルの中身を格納し、表示するプログラムです。file1の内容をfile[0].lines[○○]に、file2の内容をfile[1].lines[○○] に格納している。

出力例

出力形式の一例は以下の通り。各自、この出力例よりも優れた出力を行うこと。

Motif: MATa1

pro:YJR048W
pos:80
hit(TGATGTAATC)= 6.81

pro:YDL227C
pos:365
hit(TGATGTTAAT)= 6.51

pro:YDL227C
pos:405
hit(TGATGTAAAT)= 8.22

プログラム上のヒント

コンパイルをする手法

C言語により記述されたプログラムは直接実行できないため機械語に翻訳する作業、コンパイルが必要がある。コンパイルは以下の要領でコンソールにて行う。

gcc kadai2.c

コンパイルを行うと(a.out)という実行ファイルが同ディレクトリに配置される。このファイルの実行はfile1,file2を引数とする場合

./a.out file1 file2

といった要領で行うことができる。

読み込むファイルを引数で指定する

本課題では、プログラムを動かすごとに異なるファイル（転写因子結合部位配列ファイル）を指定することになるため、毎回プログラムを書き換えるのではなく、引数からファイル名を得る仕様にすると良い。

例えばサンプルプログラムを実行するときは、以下のように実行するプログラムprogの後ろに引数を指定する。なお、引数指定するファイルは、実行ファイルと同一のプロジェクトに入れておく。

prog filename1 filename2

ファイル名	説明
`filename1`	転写因子結合部位配列（のマルチプルアラインメント）ファイル
`filename2`	8 つのプロモータ配列ファイル

プログラム内で、引数の個数はargcに、引数のそれぞれはargv[1], argv[2], … に格納される。 argcおよびargvは以下のように記述されている（実際にサンプルプログラムの中身を見てみてください）。

#include <stdio.h>
int main(int argc, char *argv[])
{
  ....
  return 0;
}

ファイルを読み込む

ファイルに書かれている情報をプログラム内で読込むためには、以下のような関数が有用である (詳細は各自教科書を確認すること)。

関数	説明
`fgetc()`	ファイルから一文字ずつ読込みを行う
`fgets()`	ファイルから一列ずつ読込みを行う

fgets() を用いた場合のプログラム例は以下の通り。

#include <stdio.h>
#define BUFSIZE 1000
#define MAX_NAME 100
int main(int argc, char* argv[])
{
    FILE *fp;
    char buffer[BUFSIZE];
    char filename[MAX_NAME];
    ... 
    strcpy(filename, argv[1]); 
        /* コマンドラインで指定したファイル名を、
         * 「filename」という変数にコピーしている */
    ... 
    fp = fopen(filename, "r");
        /* 「filename」によって指定されているファイルを、
         * 「読込みモード」で開く */
    if(fp == NULL){
        fprintf(stderr, "%s open error?n", filename); exit(1);
    }
        /* うまくファイルを開くことができなかったとき、
         * エラー情報を出力するためのおまじない */
    while(fgets(buffer, sizeof(buffer), fp) != NULL){
        /* 一回のループごとに buffer には一行 (最大 BUFSIZE バイト)
         * の文字列が入り, 文字列末尾の改行は'?0'に変換される。 
         * ファイルの全ての列を読込むと、ループから脱出する。
         * 良く使われる構文だが、おまじないとして覚えることも可 */
        ...
        /* このループ内で何らかの処理を行うことにより、
         * ファイルから読込んだ内容を出力結果に反映させる 
         * (今回は、buffer に代入された文字列を何らかの変数に格納しておく) */ 
        ...
    }
    ...
    fclose(fp);
        /* 開いたファイルは必ず閉じること！ */
}

構造体を使う

例えば, プロモーター配列の名前や配列, モチーフ位置等は, 構造体として1つにまとめておくと便利である。

typedef struct promoter_st {
    char name[MAX_LINE];
    char seq[MAX_LINE];
} promoter;

同様に、転写因子の結合部位配列にも構造体を適用すると便利であるが、構造体のメンバにどのような要素を含めるかは各自で考える。

配列のサイズを可変にする

配列サイズは固定にするだけでなく、可変に(動的に)することも可能である。この動作を実現するためには、以下の関数が有用である (詳細は各自教科書を確認すること)。

関数	説明
`malloc()`	配列のために連続された領域を確保する
`calloc()`	配列のために非連続な領域を確保する
`free()`	`malloc()` や `calloc()` によって確保された領域を開放する

malloc() 関数を使用した場合のプログラム例は以下の通り。malloc() や calloc() によって確保された領域は何もしないとプログラムが終了するまで確保されたままなので、free()で解放することを忘れないように注意する。

void func(){
    char *array;
    int len;
    ...
    /* len に適切な値を入れる */
    ...
    array = (char *)malloc(sizeof(char) * len);
    ...
    free(array);
}

文字配列の長さを知る

C言語にはstringという概念がなく,char型の配列を扱っているため,文字列の長さを得るために関数を用いる。ある文字配列の長さを知りたい場合、strlen()という関数が有用である。strlen()を使用するために

#include <string.h>

を先頭に宣言する必要がある。

strlen()を用いたプログラム例は以下の通り。

#include <stdio.h>
#include <string.h>

main(){
  char buffer[1024];
  printf("好きなストリングを入力してください（例えばモチーフ配列）：");
  fgets(buffer,1024,stdin);buffer[strlen(buffer)-1]=0; 
        /* 入力して、最後の改行コードを消す */
  int length=strlen(buffer);
        /* 長さを図る */
  printf("入力したストリング( %s )は%d文字です。\n",buffer,length);
}

配列を関数の引数として渡す

1次元配列の場合のプログラム例を以下に示す。

#define LENGTH 10
void func(int *array){
    /* void func(int array[]) と書くこともできる */
    ....
}

int main(){
    int array[LENGTH];
    ....
    func(array);
    ....
}

2次元配列の場合のプログラム例を以下に示す。

#define ROW  5
#define COL 10
void func(int matrix[][COL]){
    /* void func(int matrix[ROW][COL]) と書くことも、
     * void func(int (*matrix)[COL]) と書くこともできる */
    ...
}

int main(){
    int matrix[ROW][COL];
    ...
    func(matrix);
    ...
}

列挙型を使用する

本プログラムでは、転写因子結合部位配列の各座位に各塩基が何回ずつ現れたかをカウントする必要がある。この場合、Aは配列の最初の要素に対応させる、といった対応付けが必要となる。上記の対応付けは、ときにプログラムのいたる所で行われ、部分部分で対応関係を間違えてしまうと、悲惨な結果を招いてしまう。そこで対応関係を一箇所に集約させる必要があるが、そのための方法論が列挙型の仕様である。列挙型を用いたプログラム例を以下に示す。

...
enum dna {A, C, G, T}; 
    /* A=0, C=1, G=2, T=3 となる */
...
void f(char *filename) 
{
    /* 塩基文字を列挙していたおかげで、
     * 'A' は「0」といった対応関係を一箇所にまとめることができる  */
    switch (buffer[i]){
        case 'A' : matrix[A][i]++; break; 
        case 'C' : matrix[C][i]++; break;
        case 'G' : matrix[G][i]++; break;
        case 'T' : matrix[T][i]++; break;
        default  : fprintf(stderr, "unknown char %c?n", buffer[i]);
    }
}

log()を使用する

今回の実験では log-odds scoreを求めるためにlog()を使用する。log()を使用するために

#include <math.h>

を先頭に宣言する必要がある。

その他のポイント

マジックナンバーを使わずに、`#define` などを使う

マジックナンバーとは 0 とか 1 とかの具体的な数字のこと。下の悪い例のように、プログラム内に直接数字を書き込むことは望ましくないので避ける。

void func(){
   int array[10];         /* 直接, 配列サイズを書き込む */
   ...
   for(i=0 ; i<10 ; i++){ /* 配列サイズをここでも使うので,
                           * 配列サイズを変更する場合, 
                           * 全てを書き直さなくてはいけない */
       ...
   }
   ...
}

#define 指定をすることで、マジックナンバーは使わずに済む。

#define LENGTH 10          /* LENGTH が 10 と同じになる
                            * 配列サイズを変更する場合, 
                            * ここだけを変更すればよい */
void func(){
   int array[LENGTH]; 
   ...
   for(i=0 ; i<LENGTH ; i++){
       ...
   }
   ...
}

デバッグをする方法

printf() どこまで実行できたかを確認するためにprintf()で適当な文字を出力する。セグメンテーションフォールトの発生位置などを知るのに有用である。
gdb セグメンテーションフォールトが起きているときにとても役に立つ。 gdbを使うときにはgcc -gでコンパイルする。とても優れたツールで何でもできるのですが、詳しい使い方は各自勉強してください。

コマンド	説明
set args `<hogehoge>`	実行する時にコマンドライン引数を `<hogehoge>`に設定する
run	プログラムを実行する
bt	backtrace. コールスタックを表示する。（呼ばれてきた関数や引数）
list	周りのソースコードを表示する
break <行/関数名>	<行/関数名>まで実行したら、一時止まる
cont	continue. 一時止ってから、また続ける
print <表現>	<表現>を計算して、表示する（例えば「print i」はiという変数の中身を表示する）
n	next. 次の行まで実行進める
s	step. 別の行が始まるまで実行進める（つまり別の関数が呼ばれたら、その関数の中に入って、また止まる）
quit または ctrl-d	gdb終了

発展課題

スコア分布の統計的な考察

補足スライドを参考に、閾値について統計学的に検討してください。

検討１：
スコアの分布を作成し、これを正規分布と仮定して検討する。
検討２：
zスコアを算出し、標準正規分布上で検討する。

プロモーター配列からコンセンサス配列を抽出する

これまではコンセンサス配列が既知であるとして実験を進めてきたが、ここではプロモーター配列からコンセンサス配列を抽出する。余裕のある人は挑戦してみてください。

ギブスサンプリングを用いる方法

補足スライドを参考に、ギブスサンプリングを用いて、転写因子REB1が結合すると知られている遺伝子群のプロモーター配列（20本）からモチーフ配列を求めよ。プロモーター配列（FASTA形式）とモチーフ長は以下の通り。

説明	ファイル名
転写因子 REB1 が結合する遺伝子群のプロモーター配列	REB1.promoters
転写因子 REB1 のモチーフ長	7

マルチプルアライメントを用いる方法

酵母の Msn2 というタンパク質は GSY2, CCP1, GPH1 等の遺伝子のプロモータ配列と結合し、それら遺伝子の転写活性を制御していることが知られている。このことからこれら3遺伝子のプロモーター配列はMsn2の結合する共通のコンセンサス配列を共有していると予想できる。そこで、三つの遺伝子のプロモーター領域のローカルマルチプルアライメントを計算し、高度に保存されている配列を発見する。

(ローカル)マルチプルアライメントを用いて、GSY2、CCP1、GPH1のプロモーター配列（3本）からMsn2の結合部位モチーフ配列を求めよ。3遺伝子の上流領域300 bpの配列（FASTA形式）は以下の通り。

説明	ファイル名
タンパク質 Msn2 が結合する遺伝子（GSY, CCP1, GPH1）のプロモーター配列	Msn2.txt

ペアワイズ・グローバルアライメント：

「バイオプログラミング第二」の講義スライドを復習せよ。

ペアワイズ・ローカルアライメント：

【表計算について】
グローバルアライメントが配列全体をアライメントするのに対し、ローカルアライメントは高度に保存された配列のみをアライメントする。そのためグローバルアライメントにおけるスコア表で、負のスコアとなってしまうようなセルには「0」を入れておき、部分配列のアライメントスコアの最大値を求められるように変更する。このことは、「バイオプログラミング第二」資料 8 ページ、(2) 式に以下の case 4 を加え、4 つの値の最大値をそのセルの値とすることで実現できる。

dp[i-1][j-1] + s(x[i], y[i]) /* case (1) */
dp[i-1][j] + d               /* case (2) */
dp[i][j-1] + d               /* case (3) */
0                            /* case (4) */

また、初期化の式は以下のように変更する。

dp[0][0] = 0
dp[i][0] = 0
dp[0][j] = 0

【トレースバックについて】
ローカルアライメントのトレースバックは、スコア表の最大の値を示すセルからスタートする。これは、グローバルアライメントではスコア表の右下からスタートすることと対照的である。

ローカル・マルチプルアライメント：

ここでは、三本の配列をマルチプルアライメントすることを考える。配列が三本あるため、ペアワイズアライメントでは二次元であったスコア表を、三次元に変更する。三次元のスコア表は、以下のように初期化される。

dp[0][0][0] = 0
dp[0][j][k] = 0
dp[i][0][k] = 0
dp[i][j][0] = 0

三次元のスコア表は、以下の計算を表の左上から繰り返すことで埋めることができる。

dp[i][j][k] = max{
 dp[i-1][j-1][k-1] + s(x[i], y[j]) + s(x[i], z[k]) + s(y[j], z[k])  /* case (1) */
 dp[i-1][j-1][k]   + s(x[i], y[j]) + d             + d              /* case (2) */
 dp[i-1][j][k-1]   + d             + s(x[i], z[k]) + d              /* case (3) */
 dp[i][j-1][k-1]   + d             + d             + s(y[j], z[k])  /* case (4) */
 dp[i-1][j][k]     + d             + d             + 0              /* case (5) */
 dp[i][j-1][k]     + d             + 0             + d              /* case (6) */
 dp[i][j][k-1]     + 0             + d             + d              /* case (7) */
 0                                                                  /* case (8) */
}

トレースバックは、ペアワイズのローカルアライメントと同様に、スコア表中の最大の値を示すセルからスタートする。強く保存されたコンセンサス配列を検出したいときには、ギャップペナルティを大きく設定すると良い。今回は、以下のパラメータを用いて、マルチプルアライメントを行うこと。

種類	スコア
マッチ	+1.0
ミスマッチ	-1.0
ギャップ	-5.0

スコア表のデータ構造：

300 × 300 × 300 の三次元配列を確保する際に、以下のように宣言したのではセグメンテーションエラーが起きてしまう可能性が高い。

int dp[300][300][300];

そのような問題が起こったならば、基礎課題のプログラムのヒントで示した方法で、動的にメモリを確保すると良い。

正解 (部分的に公開)：

「GTTT」で始まり、「GA」で終わる長さ 17 のコンセンサス配列が正解である。

さらに発展問題

発展課題で触れたようなマルチプルアライメントの方法を、動的計画法 (Dynamic Programming) と呼ぶ。動的計画法は最適解を得られるが、配列の本数が増えるにつれ、計算時間・計算メモリが指数関数的に増加してしまう。そこで配列の本数が 5 本以上の際には累進法 (Progressive Alignment) と呼ばれる方法が用いられる。例えば、ClustalW もまた、累進法によりマルチプルアライメントを計算している。累進法によりマルチプルアライメントを計算せよ。
発展課題では、アライメントする配列の本数を「三本」と決めうちにした。この決めうちをせずに、入力ファイルに含まれている配列の本数を可変とするようなプログラムを実装せよ。

レポートの書き方

「2.スコア行列を用いた転写因子結合部位の予測」については、以下の項目をレポートにまとめてください。

目的
本実験の目的を簡潔に書く。
方法
プログラムをどのようなデータ構造、アルゴリズムを用いて作成したかをわかりやすく(簡単に)書く。
結果
- プログラムの出力結果を同一プロジェクト内に置く(ファイル名は’result.txt’などとする)。教員及びTAが再現できるように入出力が分かるようにすること。
- 結合部位の予測時、スコア行列への一致度の閾値も記載する。
- 予測結果から8つのプロモータ配列と8つの転写因子間の転写制御関係図を作成する。（参考：講義スライド13ページ目）
考察
- スコア行列への一致度の閾値に対する考察（なぜ、その閾値にしたかなど）
- 予測された転写因子とプロモータ配列間の転写制御関係について、関連する論文、また関連する図書などを参考にして、生物学的な知見、また考察を述べる。
- 予測された中に偽陽性があると考えられる場合、また実際には制御関係が存在するのに正しく予測されなかったものがあると考えられる場合、なぜこれらが発生したのかを考察する。
- その他、各自オリジナルの考察 (可能であればやってください。)
結論
本実験の結論を簡潔に書く。
参考資料
参考にした論文 or 図書 or ホームページのURLを記載する。 (Wikipediaであろうと参考にしたものはすべて)

発展課題をやった人は適宜追加して書いてください。

レポート課題のヒント

出芽酵母の遺伝子には「遺伝子名」と「ORF名 (OpenReadingFrame : 染色体上の位置により機械的につけられた名前)」の2つがある。今回用意したプロモーター配列にはORF名のみ記されているが、本課題では遺伝子名に対応付けて考えた方がよい。遺伝子名を調べたり、生物学的な考察をする際に、転写因子と遺伝子の制御関係をまとめた次のようなサイトが参考になる。

SCPD (The Promoter Database of Saccharomyces cerevisiae)
SGD (Saccharomyces Genome Database)

論文はPubMedで検索できます。

概要