アミノ酸残基や糖残基が前後の残基と結合する際、取ることのできる角度の範囲を2 次元平面に図示したものをラマチャンドランプロットといいます。この分布からポリペプチド鎖や糖鎖が取りうる立体構造についての情報を考察できます。このラマチャンドランプロットを、タンパク質の構造を収録したpdbファイルを処理することにより描画するプログラムを作成しました。PDBファイルは、タンパク質と核酸の3次元構造の構造座標（立体配座）を蓄積している国際的な公共のデータベースです。PDBに蓄積されている構造データは、X線結晶解析法、NMR法（核磁気共鳴法）などによって実験的に決定されたデータです。
ラマチャンドランプロットは、ポリペプチドを考えたとき、構成単位はアミノ酸ということになります。縮合前のひとつひとつのアミノ酸にはカルボキシル基（COOH）、水酸基（OH）、水素（H）、側鎖の４つが付いている炭素があります。これをα炭素と言います。ペプチド結合することによって、・・・-α炭素-ペプチド結合-α炭素-ペプチド結合-・・・の繰り返しが出来上がります。ここで、ペプチド結合は-CONH-の４原子からなりますが、C=OとN-H間の二重結合の共鳴により平面構造を取ります。
したがって、ポリペプチドの構造はペプチド結合面どうしの二面角がすべて明らかになれば一意に定まります。定義としては、二面角は2つの交差する面の間の、交線に垂直な3つ目の面上での角度です。Cα-N結合角はφ、Cα-C結合角はψで、それぞれの角度はCαから見て時計回りに計算します。この二面角φ、ψの二次元散布図がラマチャンドランプロットです。
さて、具体的に計算する方法が書いてあるものは、検索した限りではあまりなかったのですが、平面どうしのなす角度を求めればよいので外積どうしの内積を求め、逆余弦関数で角度に戻してあげればよいでしょう。右回りか左回りかの判定は少し回りくどいのですが、外積計算で得たベクトルと平面内の回転軸となっているベクトルとの内積の正負判定で可能です。
入力としてpdbファイルを用意します。フォーマットの一部を以下に示すとおりです。ATOM定型句の次には左の列から順にN末端からの原子連番、原子の種類、その原子が構成するアミノ酸の種類、x座標、y座標、ｚ座標となっています。それより右の列は、今回は使用しません。具体的なpdbファイルとして、麦角アルカロイド生合成に関与する酵素のpdbファイル（4QNW.pdb）を用意しました。
ATOM 73 N GLY A 12 35.059 -0.373 4.881 1.00 3.78 N ATOM 74 CA GLY A 12 35.666 0.667 4.102 1.00 3.04 C ATOM 75 C GLY A 12 36.817 0.200 3.245 1.00 3.34 C ATOM 76 O GLY A 12 37.225 -0.981 3.177 1.00 2.33 O
出力はGNUPLOTと連携して散布図を描画します。同時にテキストファイルに二面角の一覧を以下の例のように出力します。
354 1.242975 -0.742149 LEU 355 -1.375875 2.770974 SER 356 -1.142880 -0.632364 ARG
角度を求めるために外積どうしの内積を計算しましたが、求まった角度が右回りか左回りかを判定するために別の計算も計算しています。ラマチャンドランプロットの描画に必要なプログラムはすでにあるのですが勉強のために作ってみました。
#include <string.h>
#include <math.h>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#define MaxResidue 500
#define NN 256
class coordinates
{
public:
double x;
double y;
double z;
};
class amino_acid
{
public:
int number;
coordinates CA;
coordinates N;
coordinates C;
char TYPE[256];
};
int number_of_residue=0;
double temp_x;
double temp_y;
double temp_z;
int aminonum;
char tAtom[256]={'\0'};
char tType[256]={'\0'};
char tChain[256]={'\0'};
int temp;
char temp_char[256]={'\0'};
amino_acid aminoacid[MaxResidue];
double phi[MaxResidue];
double psi[MaxResidue];
/* read */
int read(void) {
char ctemp_x[256]={'\0'};
char ctemp_y[256]={'\0'};
char ctemp_z[256]={'\0'};
int aminonum;
char tAtom[256]={'\0'};
char tType[256]={'\0'};
char tChain[256]={'\0'};
FILE *fp;
char *filename = "4OC6.pdb";
char readline[NN] = {'\0'};
/* ファイルのオープン */
if ((fp = fopen(filename, "r")) == NULL) {
fprintf(stderr, "%sのオープンに失敗しました.\n", filename);
exit(EXIT_FAILURE);
}
/* ファイルの終端まで文字を読み取り表示する */
while ( fgets(readline, NN, fp) != NULL ) {
if(strncmp(readline, "ATOM", 4)==0){
sscanf(readline, "%s %d %s %s %s %d %s %s %s %d %d", temp_char, &temp, tAtom, tType, tChain, &aminonum, ctemp_x,ctemp_y,ctemp_z, &temp, &temp);
aminoacid[aminonum].number=aminonum;
strcpy(aminoacid[aminonum].TYPE, tType);
if(strlen(tAtom)==2 && strncmp(tAtom, "CA", 2)==0){
number_of_residue+=1;
aminoacid[aminonum].CA.x=atof(ctemp_x);
aminoacid[aminonum].CA.y=atof(ctemp_y);
aminoacid[aminonum].CA.z=atof(ctemp_z);
} else if(strlen(tAtom)==1 && strncmp(tAtom, "N", 1)==0){
aminoacid[aminonum].N.x=atof(ctemp_x);
aminoacid[aminonum].N.y=atof(ctemp_y);
aminoacid[aminonum].N.z=atof(ctemp_z);
} else if(strlen(tAtom)==1 && strncmp(tAtom, "C", 1)==0){
aminoacid[aminonum].C.x=atof(ctemp_x);
aminoacid[aminonum].C.y=atof(ctemp_y);
aminoacid[aminonum].C.z=atof(ctemp_z);
}
}
}
fclose(fp);
return EXIT_SUCCESS;
}
coordinates vec(coordinates a,coordinates b){ //2個目に向かうベクトル
coordinates vec;
vec.x=b.x-a.x;
vec.y=b.y-a.y;
vec.z=b.z-a.z;
return vec;
}
coordinates cross_product(coordinates a,coordinates b){
coordinates vec;
vec.x = a.y*b.z - a.z*b.y;
vec.y = a.z*b.x - a.x*b.z;
vec.z = a.x*b.y - a.y*b.x;
return vec;
}
double inner_product(coordinates a,coordinates b){
double value;
value=a.x*b.x+a.y*b.y+a.z*b.z;
return value;
}
double magunitude(coordinates a){
double mag;
mag=pow(a.x*a.x+a.y*a.y+a.z*a.z,0.5);
return mag;
}
int sgn(double n){
return (n > 0) - (n < 0);;
}
/* cal phi and psi */
int phi_and_psi(int num) {
coordinates cross_product_N;
coordinates cross_product_CA;
coordinates cross_product_C;
coordinates v_PreCtoN;
coordinates v_NtoCA;
coordinates v_CAtoC;
coordinates v_CtoPostN;
v_PreCtoN=vec(aminoacid[num-1].C,aminoacid[num].N);
v_NtoCA=vec(aminoacid[num].N,aminoacid[num].CA);
v_CAtoC=vec(aminoacid[num].CA,aminoacid[num].C);
v_CtoPostN=vec(aminoacid[num].C,aminoacid[num+1].N);
cross_product_N=cross_product(v_PreCtoN,v_NtoCA);
cross_product_CA=cross_product(v_NtoCA,v_CAtoC);
cross_product_C=cross_product(v_CAtoC,v_CtoPostN);
phi[num]=sgn(inner_product(v_PreCtoN,cross_product_CA))*acos(inner_product(cross_product_N,cross_product_CA)/(magunitude(cross_product_N)*magunitude(cross_product_CA)));
psi[num]=sgn(inner_product(v_NtoCA,cross_product_C))*acos(inner_product(cross_product_CA,cross_product_C)/(magunitude(cross_product_CA)*magunitude(cross_product_C)));
std::cout<<num<<" "<<phi[num]<<" "<<psi[num]<<std::endl;
return 0;
}
void exportData(void){
FILE *fp ;
if((fp=fopen("/Users/YM/Desktop/save.txt","w"))!=NULL){
for(int i=1;i<number_of_residue;i++){
fprintf(fp,"%d %f %f ",i,phi[i],psi[i]);
fputs(aminoacid[i].TYPE, fp);
fprintf(fp,"\n");
}
fclose(fp);
}
}
int plot()//パイプによってグラフ描画
{
FILE *gp;
/* ---- パイプを開く ---- */
gp = popen("/usr/local/bin/gnuplot -persist","w");
fprintf(gp, "set terminal x11\n");
/* ---- グラフ作成 ---- */
fprintf(gp, "set xtics 1.5707\n");
fprintf(gp, "set ytics 1.5707\n");
fprintf(gp, "set size ratio 1\n");
fprintf(gp, "set xlabel \"phi\"\n");
fprintf(gp, "set ylabel \"psi\"\n");
fprintf(gp, "plot \"save.txt\" using 2:3\n");
fprintf(gp,"\n");
pclose(gp);
return 0;
}
/* main */
int main(void) {
read();
std::cout<<number_of_residue<<std::endl;
for(int i=2;i<number_of_residue-1;i++){
phi_and_psi(i);
std::cout<<i<<" "<<phi[i]<<" "<<psi[i]<<std::endl;
}
exportData();
plot();
}