接触原子の計算
今回はタンパク質の立体構造座標を使った計算を行うプログラムを作ります。具体的には、リガンドとの結合に関わるタンパク質原子を求めるプログラムです。タンパク質を構成する原子とリガンドを構成する原子の距離が4.0Å未満（水分子が入り込めない距離）の場合にそれらの原子は接触していると見なし、接触している原子ペアをリストアップしてみましょう。
このプログラムで必要なことは、
- PDBの座標データからタンパク質の原子座標データを読み込む。
- PDBの座標データからリガンドの原子座標データを読み込む。
- タンパク質原子とリガンド原子の距離を計算する。
- 原子間距離が4.0Å未満の原子ペアを出力する。
の４つの行程です。
このうち、原子座標のデータですが、純粋に距離だけを計算するだけであればxyz座標だけを読み込めば十分です。しかし、今回作りたいプログラムでは最後にどの原子ペアが接触しているかを出力するので、xyz座標がどの原子のものなのかも記録しておく必要があります。一番単純にはxyz座標と原子名、アミノ酸残基名などを記憶する変数を別々に作成する方法がありますが、それではいろいろと管理が面倒です。今回の例のようなxyz座標と原子名などの異質なデータを一度に管理する仕組みがＣ言語にはあり、構造体と呼ばれます。構造体は、
struct HOGE {
int member1;
char member2[5];
float member3;
} hoge;
という形で宣言されることが基本となります。この宣言の場合、メンバーとしてint型の変数member1とchar型の配列member2[5]、float型の変数member3が含まれているHOGEという構造体としてhogeが作られます。それぞれのメンバーにアクセスする場合は、例えば、
printf("%d %s %f\n", hoge.member1, hoge.member2, hoge.member3);
のようになり、構造体の名前とメンバーの名前を”.”（ピリオド）でつなぎます。
なお、構造体の宣言でよく使われている方法として、以下のようなものがあります。
typedef struct HOGE {
int member1;
char member2[5];
float member3;
} hoge_t;
詳しい説明は省きますが、このように宣言することで”hoge_t”という文字列にintやcharと同じような役目を持たせることができ、この宣言の後に
hoge_t hoge;
とさらに宣言すると、HOGE構造体であるhogeが宣言できます。いってしまえば、hoge_t型という変数ができたようなものです。
では、この構造体の話をふまえてプログラムの中身を考えてみましょう。
main()関数
まず、各種の準備とmain()関数の中身は以下のようになります。このプログラムでは、コマンドラインの引数からPDBファイルの名前を受け取るように設計しています。
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#define MAXCHR 128
#define MAXATOM 5000
#define MAXHET 200
typedef struct PDB {
char atomName[5];
char resName[4];
char ch;
char resNum[5];
float e[3];
} pdb_t;
int main(int argc, char *argv[])
{
int atomNum;
int hetNum;
pdb_t atom[MAXATOM];
pdb_t hetatm[MAXHET];
float dis[MAXATOM][MAXHET];
FILE *fp;
int ReadPDB(FILE *, pdb_t [], char *);
void CalcDis(int, int, pdb_t [], pdb_t [], float [][MAXHET]);
void ShowContact(int, int, pdb_t [], pdb_t [], float [][MAXHET]);
if (argc != 2) {
fprintf(stderr, "%s PDBFILE\n", argv[0]);
exit(1);
}
if ((fp = fopen( argv[1], "r")) == NULL) {
fprintf(stderr, "main(): Cannot open file %s\n", argv[1]);
exit(1);
}
atomNum = ReadPDB(fp, atom, "ATOM ");
rewind(fp);
hetNum = ReadPDB(fp, hetatm, "HETATM");
CalcDis(atomNum, hetNum, atom, hetatm, dis);
ShowContact(atomNum, hetNum, atom, hetatm, dis);
return 0;
}
10行目から16行目に先ほど説明した構造体の宣言があります。この宣言をこの場所ですることで、これ以降のソースコードでは”pdb_t”という文字列があれば、atomName[5], resName[4], ch, resNum[5], float[3]のメンバーを持つ構造体型を示します。一般には、typedefによる構造体の宣言は、関数の外で、ソースコードの最初にまとめて行います。
22行目と23行目でpdb_t型の構造体を宣言しています。それぞれ、タンパク質を構成する原子（以降、atom原子とします）のデータをしまうpdb[]とリガンドを構成する原子（以降、hetatm原子とします）のデータをしまう構造体です。
27行目から29行目で使う関数の宣言をしていますが少し注意点があります。CalcDis()関数およびShowContact()関数に渡す引数の最後に”float [][MAXHET]”という記述があります。これは、２次元の配列を引数として渡すことを宣言する場合は、最後の次元の大きさを必ず指定する必要があるためです。
atom原子とhetatm原子のデータは41行目と43行目のReadPDB()関数で読み込みます。PDBのデータを思い出してもらうとわかりますが、タンパク質の原子座標の文頭は”ATOM “であり、リガンドの原子座標の文頭は”HETATM”である以外は基本的に同じフォーマットで書かれています。そのため、”ATOM “から始まる行を読むか”HETATM”から始まる行を読むかをReadPDB()関数に渡してやることで、atom原子とhetatm原子のどちらのデータでも同じ関数で読み込めるようにしています。
新しい要素として、42行目のrewind()関数があります。この関数は最後まで読み込み終えたファイルをまた最初から読み込むようにする関数です。これにより、ReadPDB()関数でatom原子を読む際に最後まで読み終えてしまったPDBのデータが入ったファイルを、もう一度最初から読めるようになります。
ReadPDB()関数
次は、PDBの座標データを読み込むための関数、ReadPDB()関数です。
ReadPDB()関数は、PDBファイルのポインタ、座標データを入れる構造体データ、どの行を読み込むかを判断するタグの３つを引数とします。返り値は、読み込んだ原子の数です。
int ReadPDB(FILE *fp, pdb_t pdb[], char *str)
{
int num = 0;
char line[MAXCHR];
while (fgets(line, MAXCHR, fp) != NULL) {
if (strncmp(line, str, 6) == 0) {
strncpy(pdb[num].atomName, &line[12], 4);
pdb[num].atomName[4] = '\0';
strncpy(pdb[num].resName, &line[17], 3);
pdb[num].resName[3] = '\0';
pdb[num].ch = line[21];
strncpy(pdb[num].resNum, &line[22], 4);
pdb[num].resNum[4] = '\0';
pdb[num].e[0] = atof(&line[30]);
pdb[num].e[1] = atof(&line[38]);
pdb[num].e[2] = atof(&line[46]);
num++;
}
}
return num;
}
6行目はこれまでにも出てきたように、ファイルfpから１行読み込み、MAXCHRの数だけ文字列lineに代入します。
7行目では文字列lineの先頭から６文字とstrの中身を比較します。strは”ATOM “もしくは”HETATM”なので、”ATOM “の場合はPDBのatom原子の座標データが書いてある行とマッチし、”HETATM”の場合はPDBのhetatm原子の座標データが書いてある行とマッチします。
8行目から19行目はごちゃごちゃと書かれていますが、pdb_t型の構造体に座標データを代入している部分です。8行目は、文字列lineの13文字目から4文字分をnum番目の構造体pdbのatomNameという文字列メンバーにコピーします。9行目はatomNameに代入した文字列の最後を明示するための命令です。この行は忘れやすいので注意して下さい。12行目を除く10行目から14行目は同じような内容です。15行目は文字列lineの31文字目から存在する数字（これはx座標になります）をfloat型の変数に変換してnum番目の構造体pdbのe[0]メンバーに代入しています。
CalcDis()関数
次は、atom原子とhetatm原子のすべてのペアについて原子間距離を計算するための関数、CalcDis()関数です。
CalcDis()関数は、atom原子の数、hetatm原子の数、atom原子の構造体データ、hetatm原子の構造体データ、原子間距離を入れる二次元配列を引数として受け取ります。この関数では特に返り値はもうけません。
void CalcDis(int atomNum, int hetNum, pdb_t atom[], pdb_t het[], float dis[][MAXHET])
{
int i, j;
float Distance(float [], float []);
for (i = 0; i < atomNum; i++) {
for (j = 0; j < hetNum; j++) {
dis[i][j] = Distance(atom[i].e, het[j].e);
}
}
}
5行目で実際の距離計算を行う関数Distance()の宣言をしています。
7行目から11行目で、すべてのatom原子とすべてのhetatm原子の間の距離を、Distance()関数にそれぞれの座標データを渡すことで計算しています。
Distance()関数
次は、３次元座標から距離を計算するための関数、Distance()関数です。
Distance()関数は、引数として３つの要素からなるfloat配列を２個受け取ります。３つの要素は、それぞれ原子のxyz座標が入ります。この関数が返す返り値は計算した距離になります。
float Distance(float e1[3], float e2[3])
{
float dis;
dis = pow(e1[0]-e2[0], 2.0) +
pow(e1[1]-e2[1], 2.0) +
pow(e1[2]-e2[2], 2.0);
return sqrt(dis);
}
(x1, y1, z1)と(x2, y2, z2)の距離dは(x1-x2)2+(y1-y2)2+(z1-z2)2の平方根を取ったものなので、まずは５行目から７行目にあるような式を計算します。pow()関数は、double型の引数を２つ取り、１つ目の引数を２つ目の引数乗する関数です。なお、pow()関数を使うにはmath.hが必要です。
ShowContact()関数
次は、原子間距離のデータから接触原子ペアを出力するための関数、ShowContact()関数です。
ShowContact()関数は、atom原子の数、hetatm原子の数、atom原子の構造体データ、hetatm原子の構造体データ、原子間距離が入った二次元配列を引数として受け取ります。
void ShowContact(int atomNum, int hetNum, pdb_t atom[], pdb_t het[], float dis[][MAXHET])
{
int i, j;
for (i = 0; i < atomNum; i++) {
for (j = 0; j < hetNum; j++) {
if (dis[i][j] < 4.0) {
printf("%3s %c %4s %5s - %3s %c %4s %5s %6.2f\n",
atom[i].resName, atom[i].ch,
atom[i].resNum, atom[i].atomName,
het[j].resName, het[j].ch,
het[j].resNum, het[j].atomName,
dis[i][j]);
}
}
}
}
この関数では特に新しいことはやっていません。i番目のatom原子とj番目のhetatm原子の原子間距離が4.0Åであれば、i番目のatom原子のデータ（原子名や残基名など）とj番目のhetatm原子のデータを出力しています。
これまでに比べて大分重たいソースコードになりましたが、以上で解説は終わりです。ソースコードにcontact.cと名前を付けて保存しましょう。次に、
cc -o contact contact.c -lm
とコンパイルします。pow()関数を使っているので-lmオプションをつけて下さい。
次にPDBから適当なファイル（例えばPDB IDが1mbs）をダウンロードして、プログラムを実行してみて下さい。うまく結果が出力されたでしょうか？segmentation faultなどになったら、ソースコードをもう一度見直してみましょう。
なお、今回作成したプログラムは、すべてのPDBデータに対応できるようにはなっていません。例えば、atom原子の数が5000を超えるタンパク質の計算しようとするとsegmentation faultになってしまうでしょう。これは、必要なメモリを割り当てる関数であるmalloc()関数を使うと解決できます。
練習問題
- 今回作成したプログラムは、水の原子と接触するタンパク質原子もリストアップする。しかし、水の原子データはPDBファイル中に大量に入っている場合があり、結果が煩雑となってしまう。プログラムを改変して、水の原子との接触は出力しないようにしなさい。
- ２つのサブユニットの間で接触している原子を求めるプログラムを作成しなさい。