C言語講座　第９回

hydropathy plot

今回作成するプログラムは、膜タンパク質の膜貫通部位を予測するためのhydropathy plotです。

膜タンパク質において膜に埋まっている部位は、その疎水的な環境に合わせて疎水性の高いアミノ酸残基が長く連続して存在しています。Kyte & Doolittleによりhydropathy indexが作成され（Kyte J & Doolittle RF, 1982, J Mol Biol, 157:105-132.）、各アミノ酸の疎水性が定量的に示されています。なお、hydropathy indexは値が大きくなるほど疎水的なアミノ酸を示し、具体的な数値は以下の表の通りです。

residue	index	residue	index
Ile	4.5	Trp	-0.9
Val	4.2	Tyr	-1.3
Leu	3.8	Pro	-1.6
Phe	2.8	His	-3.2
Cys	2.5	Glu	-3.5
Met	1.9	Gln	-3.5
Ala	1.8	Asp	-3.5
Gly	-0.4	Asn	-3.5
Thr	-0.7	Lys	-3.9
Ser	-0.8	Arg	-4.5

hydropathy indexを用いて、膜貫通部位は以下のような考え方で求めます。この考え方の肝は、１つのアミノ酸残基だけの疎水性を見るのではなく、前後のアミノ酸残基の疎水性も考慮に入れるところです。

1. 各アミノ酸残基にhydropathy indexを割り振る。
2. i番目のアミノ酸残基の前後３残基（i-3からi+3まで）のhydropathy indexの平均を求め、その値をi番目のアミノ酸残基の疎水性とする。
3. 各アミノ酸残基の疎水性が-0.4以上の場合、その残基は疎水的な領域を作る候補とする。
4. 膜貫通ヘリックスを形成できる残基数以上に疎水的な領域が存在する場合は、その部分が膜貫通部位であると考えられる。

以上の考え方に従って、FASTA形式のアミノ酸配列データを読み込んで、hydropathy plotを表示するプログラムを作成してみましょう。

まずは、FASTA形式を読み込む関数を作りましょう。この関数に渡す引数は、FASTA形式のデータが入ったファイル名とアミノ酸配列を格納するchar型の配列です。この関数の返り値は、読み込んだアミノ酸配列の長さとすると便利です。

int  ReadFASTA(char *file, char seq[])
{
    int    len = 0;
    int    flag = 0;
    char   amino;
    FILE  *fp;

    if ((fp = fopen(file, "r")) == NULL) {
        fprintf(stderr, "ReadFASTA(): Cannot find %s\n", file);
        exit(1);
    }

    while ((amino = fgetc(fp)) != EOF) {
        if (amino == '>') {
            flag = 1;
        } else if (flag == 1 && amino == '\n') {
            flag = 0;
        }

        if (flag == 0 && amino != '\n' && amino != ' ') {
            seq[len++] = amino;
        }
    }
    seq[len] = '\0';

    return len;
}

8行目から11行目でFASTA形式のデータが入ったファイルを開いています。ここで少し新しい部分は、もし指定されたファイルが存在しない場合、エラー出力をしてプログラムを終了させるところです。fopen()関数は、ファイルが開けない場合”NULL”を返します。そのため、8行目においてFILEポインタであるfpにNULLが代入されたら、9行目と10行目を実行します。9行目のfprintf()関数は、使い方はprintf()関数とほぼ同じですが、１つ目の引数にどこに出力するかを指定します。8行目の場合はエラー出力をしたいので、”stderr”を指定しています。

13行目から24行目がファイルからFASTA形式を読み込む本体となります。

この部分で使われている変数の説明をします。aminoは、ファイルに納められている文字を１つをいれる変数です。lenは、読み込んだアミノ酸配列の長さを示します。flagは、FASTA形式のヘッダー行（>で始まる行）を読み込んでいるところかどうかを判定する変数で、1だったらヘッダー行を読み込んでおり、0だったらそれ以外を読み込んでいることを示します。

13行目にあるfgetc()関数は、ファイルから文字を１つだけ読み込みます。もしファイルの最後に来たらfgetc()関数は”EOF”を返します。よって、13行目を読み下すと、「ファイルの中身を一文字ずつ読み込んでaminoに入れることを、ファイルの最後まで繰り返しなさい」となります。

14行目から18行目は、ヘッダー行かどうかを判定してflag変数を1に変えるところです。

20行目から22行目は、読み込んだ文字がヘッダー行のものではなく、かつ、改行や空白でない場合に、配列seq[]のlen番目にaminoを代入します。代入し終えたらlenに1を足しています。

24行目では、文字列の最後に文字列の終わりの印である’\0’を代入しています。

26行目でこの関数の返り値である読み込んだアミノ酸配列の長さを返しています。

次は、読み込んだアミノ酸配列のアミノ酸残基それぞれにhydropathy indexを割り振る関数を作ります。考え方の1.にあたる関数です。この関数の引数は、アミノ酸配列の長さ、アミノ酸配列が入った配列、割り振るhydropathy indexを入れる配列とします。この関数から特に返り値は必要としないので、void型の関数としましょう。void型の関数とは、返り値を返さない関数のことです。

void  AssignHydropathyIndex(int len, char seq[], float index[])
{
    int   i, j;
    char  amino[] = "IVLFCMAGTSWYPHEQDNKR";
    float HI[] = {  4.5,  4.2,  3.8,  2.8,  2.5,
                    1.9,  1.8, -0.4, -0.7, -0.8,
                   -0.9, -1.3, -1.6, -3.2, -3.5,
                   -3.5, -3.5, -3.5, -3.9, -4.5 };

    for (i = 0; i < len; i++) {
        for (j = 0; j < 20; j++) {
            if (seq[i] == amino[j]) {
                index[i] = HI[j];
                break;
            }
        }
    }
}

4行目では20種類のアミノ酸でamino[]という配列を初期化しています。5行目から8行目では、4行目で初期化したアミノ酸の順番で、HI[]というfloat配列をhydropathy indexにより初期化しています。

10から17行目は、これまでにも何度か出てきた形です。ここの部分で、各アミノ酸に対応するhydropathy indexをindex[]というfloat配列に入れています。

次は、考え方の2.にあたる関数を作ります。つまり、i番目のアミノ酸残基の前後３残基（i-3からi+3まで）のhydropathy indexの平均を求め、その値をi番目のアミノ酸残基の疎水性とする関数です。この関数の引数は、アミノ酸配列の長さ、割り当てられたhydropathy indexが入ったfloat配列、hydropathy indexの平均を入れるfloat配列となります。この関数では特に値を返さないため、void関数になっています。

void  CalcHydropathy(int len, float index[], float hydropathy[])
{
    int   i, j;

    for (i = 3; i < len-3; i++) {
        hydropathy[i] = 0.0;
        for (j = -3; j <= 3; j++) {
            hydropathy[i] += index[i+j];
        }
        hydropathy[i] /= 7.0;
    }
}

2.の考え方でちょっとした問題となるのは、前後３残基のhydropathy indexが必要となるためにN末端の３残基とC末端の３残基には平均hydropathy indexを求められないことです。そこでこのプログラムでは、N末端の３残基とC末端の３残基は計算しないようにしています。また、平均値を計算するのに、いちいちi-3番目、i-2番目、…、i+2番目、i+3番目のアミノ酸残基のhydropathy indexを書き出して合計するのではなく、for文を使って楽をしています。

次の関数は、結果の出力をする関数です。この関数の引数は、アミノ酸配列の長さ、アミノ酸配列を入れたchar配列、平均hydropathy indexを入れたfloat配列です。この関数も特に値を返さないのでvoid関数です。

void PrintHydropathy(int len, char seq[], float hydropathy[])
{
    int   i;

    for (i = 0; i < len; i++) {
        printf("%3d  %c  %5.1f  ", i+1, seq[i], hydropathy[i]);
        if (hydropathy[i] >= -0.4) {
            printf("*");
        } else {
            printf(" ");
        }
        printf("\n");
    }
}

この関数では特に目新しいことはしていません。考え方3.に習って、平均hydropathy indexが-0.4以上の場合に印（この場合は*）を表示するようにしているだけです。

最後は、これまでの関数を呼び出すmain()関数です。コマンドラインからの引数でFASTA形式のアミノ酸配列が入ったファイル名を指定するようにしましょう。

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

#define MAXLEN      2000

int main(int argc, char *argv[]) 
{
    int   seqLen;
    char  seq[MAXLEN];
    float index[MAXLEN];
    float hydropathy[MAXLEN];

    int   ReadFASTA(char *, char []);
    void  AssignHydropathyIndex(int, char [], float []);
    void  CalcHydropathy(int, float [], float []);
    void  PrintHydropathy(int, char [], float []);
    
    if (argc != 2) {
        fprintf(stderr, "%s  FASTA_FILE\n", argv[0]);
        exit(1);
    }

    seqLen = ReadFASTA(argv[1], seq);
    AssignHydropathyIndex(seqLen, seq, index);
    CalcHydropathy(seqLen, index, hydropathy);
    PrintHydropathy(seqLen, seq, hydropathy);

    return 0;
}

1行目から5行目まではheaderファイルの指定をしているところです。この場合は便宜上main()関数と一緒に書いていますが、実際のソースコードではファイルの先頭に書くようにして下さい。

9行目から12行目で、変数を宣言しています。MAXLENは#defineにより2000に置き換えられます。つまり、このプログラムでは読み込めるアミノ酸配列の最大残基長が2000未満ということを意味します。

14行目から17行目は関数の宣言をしています。宣言文では、引数がどういう変数かがわかればいいので、変数名は省略することができます。

19行目から22行目は、引数を１つ指定しないと使い方を表示して終了させる部分です。24行目から27行目は順番に関数を呼び出しています。

以上で、最低限必要なソースコードはすべて用意することができました。すべてを入力し、hydro.cとでも名前を付けてファイルに保存し、コンパイルしましょう。コンパイルできたら、実際に適当なアミノ酸配列（例えば、Swiss-ProtエントリーのOPSD_HUMANやHBA_HUMANなどをFASTA形式で保存したデータなど）を用意してプログラムを実行してみましょう。

練習問題

1. 考え方の4.に従って、膜貫通ヘリックスと予測されるアミノ酸残基に印を付けるようにしなさい。ここで、脂質二重膜を貫通するのに必要な最低アミノ酸残基長は20残基とする。
2. Swiss-Prot形式のアミノ酸配列データを読み込んでhydropathy plotを求めるプログラムを作成しなさい。