クリスタルコンタクトを求めるプログラム「symia2csv」

鬱症状が酷くて実家で療養中、研究をまとめる気も起こらず院を中退しようかなと考えている今日この頃だが、一応研究に関係する話題。

解いたタンパク質の結晶構造について、先輩から「クリスタルコンタクトを考慮した方が良いが、お前関係する残基をリストアップするプログラム書けるんじゃね」と言われたので、まず Ruby で書いてみた。しかし遅すぎたので改めて C で書いてみた。今度は十分速くなったので、公開してみます。

使い方は、結晶構造の PDB と対称分子のみを含む PDB（PyMOL などで作る）を用意して

./symia2csv PDB 対称分子のPDB 相互作用と見なす距離（Å単位）

でいけます。リダイレクトで CSV ファイルに出力すると見やすいです。 ビルドは

gcc symia2csv.c -lm -o symia2csv

あたりで。

ダウンロード：http://www.li-sa.jp/ocha3/c/symia2csv.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
 
#define LINE_MAX 100
 
typedef struct {
    char chain[2];
    char residue[4];
    unsigned short residue_num;
    char name[5];
    float x;
    float y;
    float z;
} atom_t;
 
typedef struct atom_list {
    atom_t *atom;
    struct atom_list *next;
} atom_list_t;
 
typedef struct atoms_distance {
    atom_t *a1;
    atom_t *a2;
    float distance;
 
    struct atoms_distance *next;
} atoms_distance_t;
 
atom_t* create_atom(char *line);
atom_list_t* create_atom_list(atom_t *a);
void release_atom_list(atom_list_t *start);
atoms_distance_t* create_atoms_distance(atom_t *a1, atom_t *a2, float distance);
void release_atoms_distance(atoms_distance_t *start);
float calculate_distance(atom_t *a1, atom_t *a2);
void print_atoms_distance(atoms_distance_t *ad);
 
float max_distance;
int main(int argc, char *argv[]) {
    if (argc != 4) {
        fputs("Usage: symia2csv PDB PDB_SYM MAX_DISTANCE\n", stderr);
        return 1;
    }
 
    FILE *fp_pdb, *fp_sym;
 
    if ((fp_pdb = fopen(argv[1], "r")) == NULL) {
        fprintf(stderr, "Can't read PDB file: %s\n", argv[1]);
        return 1;
    }
 
    if ((fp_sym = fopen(argv[2], "r")) == NULL) {
        fprintf(stderr, "Can't read PDB file (symmetric): %s\n", argv[2]);
        return 1;
    }
 
    if (!((max_distance = atof(argv[3])) > 0)) {
        fputs("MAX_DISTANCE must be greater than 0.\n", stderr);
        return 1;
    }
 
    char line[LINE_MAX];
    char header[7];
    unsigned pdb_atoms = 0;
    atom_list_t *atom_list_pdb_start = NULL;
    atom_t *atom;
    atom_list_t *atom_list = NULL;
    atom_list_t *next = NULL;
    for (; fgets(line, LINE_MAX, fp_pdb); atom_list = next) {
        sscanf(line, "%6s", header);
        if (!strcmp(header, "ATOM")) {
            pdb_atoms++;
 
            atom = create_atom(line);
            next = create_atom_list(atom);
 
            if (atom_list_pdb_start == NULL)
                atom_list_pdb_start = next;
            if (atom_list != NULL)
                atom_list->next = next;
        }
    }
 
    atom_list_t *atom_list_sym_start = NULL;
    atom_list = NULL;
    next = NULL;
    for (; fgets(line, LINE_MAX, fp_sym); atom_list = next) {
        sscanf(line, "%6s", header);
        if (!strcmp(header, "ATOM")) {
            atom = create_atom(line);
            next = create_atom_list(atom);
 
            if (atom_list_sym_start == NULL)
                atom_list_sym_start = next;
            if (atom_list != NULL)
                atom_list->next = next;
        }
    }
 
    fclose(fp_pdb);
    fclose(fp_sym);
 
    atom_list_t *al1, *al2;
    atoms_distance_t *ad_start = NULL;
    atoms_distance_t *ad = NULL;
    atoms_distance_t *ad_next = NULL;
    unsigned count;
    float distance;
    for (al1 = atom_list_pdb_start, count = 1; al1; al1 = al1->next, count++) {
        fprintf(stderr, "\rCalculating distances to atom %u/%u",
                count, pdb_atoms);
 
        for (al2 = atom_list_sym_start; al2; al2 = al2->next) {
            distance = calculate_distance(al1->atom, al2->atom);
            if (distance <= max_distance) {
                ad_next = create_atoms_distance(al1->atom, al2->atom, distance);
 
                if (ad_start == NULL) {
                    ad_start = ad_next;
                }
                if (ad != NULL)
                    ad->next = ad_next;
                ad = ad_next;
            }
        }
    }
    fputs("\n", stderr);
 
    printf("Chain,Residue,Residue.Num,Atom,Chain,Residue,Residue.Num,Atom,Distance\n");
    if (ad_start != NULL) {
        for (ad = ad_start; ad; ad = ad->next)
            print_atoms_distance(ad);
    }
 
    release_atom_list(atom_list_pdb_start);
    release_atom_list(atom_list_sym_start);
    release_atoms_distance(ad_start);
 
    return 0;
}
 
atom_t* create_atom(char *line) {
    atom_t *a = malloc(sizeof(atom_t));
 
    sscanf(line + 21, "%1s", a->chain);
    sscanf(line + 17, "%3s", a->residue);
    sscanf(line + 22, "%hu", &a->residue_num);
    sscanf(line + 12, "%4s", a->name);
    sscanf(line + 30, "%f", &a->x);
    sscanf(line + 38, "%f", &a->y);
    sscanf(line + 46, "%f", &a->z);
 
    return a;
}
 
atom_list_t* create_atom_list(atom_t *a) {
    atom_list_t *al = malloc(sizeof(atom_list_t));
    al->atom = a;
    al->next = NULL;
 
    return al;
}
 
void release_atom_list(atom_list_t *start) {
    atom_list_t *al = start;
    atom_list_t *next = NULL;
 
    for (; al; al = next) {
        next = al->next;
        free(al->atom);
        free(al);
    }
}
 
atoms_distance_t* create_atoms_distance(atom_t *a1, atom_t *a2, float distance) {
    atoms_distance_t *ad = malloc(sizeof(atoms_distance_t));
    ad->a1 = a1;
    ad->a2 = a2;
    ad->distance = distance;
 
    ad->next = NULL;
 
    return ad;
}
 
void release_atoms_distance(atoms_distance_t *start) {
    atoms_distance_t *ad = start;
    atoms_distance_t *next = NULL;
 
    for (; ad; ad = next) {
        next = ad->next;
        free(ad);
    }
}
 
float calculate_distance(atom_t *a1, atom_t *a2) {
    float dx = a2->x - a1->x;
    float dy = a2->y - a1->y;
    float dz = a2->z - a1->z;
 
    return sqrt(dx * dx + dy * dy + dz * dz);
}
 
void print_atoms_distance(atoms_distance_t *ad) {
    printf("%s,%s,%hu,%s,",
            ad->a1->chain,
            ad->a1->residue,
            ad->a1->residue_num,
            ad->a1->name);
    printf("%s,%s,%hu,%s,",
            ad->a2->chain,
            ad->a2->residue,
            ad->a2->residue_num,
            ad->a2->name);
    printf("%f\n", ad->distance);
}

malloc で確保した領域のチェックを行っていないのは良くないが、よほどメモリが足りなくなければ動くと思うのでちょっと放置。研究の役に立てば幸いです。