2010年01月04日

DNA鎖をコドン変換する簡単コマンド

仕事とあまり関係ないのですが、バクテリア等の実験で得られたDNA鎖をコドン表という表に従ってペプチドの配列に変換したいことがあります。
多分よくやられる変換だと思うのですが、そのようなプログラムが
google検索ですぐに見つからなかったので自分で書いてみました。
中身はあまり洗練されていませんが、よかったらお使いください。
使用、再配布自由です。

gcc -o dna2codon dna2codon.c

でコンパイルし、

./dna2codon ATCGGGCTTACC

でペプチド配列が表示されます。
また、途中にストップコドンがある場合はそこでカットするようにしています。
大量の配列を変換する場合は、シェルスクリプトの中でこれを使ってください。

ちなみにコドン表はここの表を使いました。
RNAバージョンもここのテーブルを使って同じように作れると思います。

///////////////////////////////////////////////////////////////////
// @create 2009.12.6
// @author t.mimori @ RAXUS Co. Ltd.
// @name dna2codon
// @stdin DNA鎖
//
// @stdout codon配列
// @brief
// DNA鎖をcodon配列に変換するプログラム
// MAX_CODON × 3個以内に限る
//
///////////////////////////////////////////////////////////////////

#define MAX_CODON 2000
//#define DEBUG

#include <stdio.h>
#include <ctype.h> //tolower()で必要

typedef enum {inv = -1, t = 0, c, a, g} tcag;

tcag baseNum(char base)
{
switch ((int) base)
{
case 't':
return t;
case 'c':
return c;
case 'a':
return a;
case 'g':
return g;
default:
return inv;
}
return inv;
}

int getCodonNum(const char *bases)
{
//計算の順番に注意
return baseNum(bases[2])
+ baseNum(bases[1]) * 4 * 4
+ baseNum(bases[0]) * 4;
}

char getCodonNameFromNum(int codonNum)
{
if (codonNum <= getCodonNum("ttc")) return 'F';
if (codonNum <= getCodonNum("ctg")) return 'L';
if (codonNum <= getCodonNum("ata")) return 'I';
if (codonNum <= getCodonNum("atg")) return 'M'; //開始
if (codonNum <= getCodonNum("gtg")) return 'V';
if (codonNum <= getCodonNum("tcg")) return 'S';
if (codonNum <= getCodonNum("ccg")) return 'P';
if (codonNum <= getCodonNum("acg")) return 'T';
if (codonNum <= getCodonNum("gcg")) return 'A';
if (codonNum <= getCodonNum("tac")) return 'Y';
if (codonNum <= getCodonNum("tag")) return '$'; //終了codon(出力されない)
if (codonNum <= getCodonNum("cac")) return 'H';
if (codonNum <= getCodonNum("cag")) return 'Q';
if (codonNum <= getCodonNum("aac")) return 'N';
if (codonNum <= getCodonNum("aag")) return 'K';
if (codonNum <= getCodonNum("gac")) return 'D';
if (codonNum <= getCodonNum("gag")) return 'E';
if (codonNum <= getCodonNum("tgc")) return 'C';
if (codonNum <= getCodonNum("tga")) return '$';
if (codonNum <= getCodonNum("tgg")) return 'W';
if (codonNum <= getCodonNum("cgg")) return 'R';
if (codonNum <= getCodonNum("agc")) return 'S';
if (codonNum <= getCodonNum("agg")) return 'R';
if (codonNum <= getCodonNum("ggg")) return 'G';
return '-'; //Error
}

int output(const char *codons)
{
printf("%s", codons);
return 0;
}

int main(int argc, char *argv[])
{
//引数の数をチェック
if (argc != 2)
{
return -1;
}

char *dna = argv[1];

int i, j, k;
i = j = k = 0;
char codons[MAX_CODON] = ""; //出力codon配列
char unitCodon[4] = ""; //1codonとなる3つの塩基

while (dna[i])
{

#ifdef DEBUG
printf("dna[%d]: %c\n", i, dna[i]);
#endif
char base = tolower(dna[i]);
if (baseNum(base) < 0)
{
return output(codons);
}
else
{
unitCodon[j] = base;
j++;
}

if (j == 3)
{
#ifdef DEBUG
printf("unitCodon: %s\n", unitCodon);
#endif
char codon = getCodonNameFromNum(getCodonNum(unitCodon));
if (codon == '$') return output(codons);
else codons[k] = codon;
#ifdef DEBUG
printf("codons: %s\n", codons);
#endif
k++;
j = 0;
}

i++;


return output(codons);
}


タグ:leopard
posted by ラクサス at 09:32| Comment(1) | TrackBack(0) | 【プログラミング】C編 | このブログの読者になる | 更新情報をチェックする
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。