情報基礎A 「Cプログラミング」総合演習
このページでは、塩基配列データをCのプログラムで扱うための具体的な手順を説明する。
プログラミングのための準備
塩基配列データの入手
現在までに知られている塩基配列の多くは、データベースを使って簡単に検索することができる。 例えば、 Genbank(NCBI)のヌクレオチド検索データベース、 日本DNAデータバンク(DDBJ)の検索ページ、 EMBL-EBIのデータベース、 などは良く知られている。 それらの中から、必要なデータを絞り込むには、専門的な知識や用語を知っておく必要があるので、 こうしたデータベースを自在に使いこなせるようになには、専門の勉強を積む必要があるだろう (何せ、この講義の担当教員の専門は物理なので、「教養」程度のことしか分からない)。
そうは言っても、検索自体は、簡単に行うことができる。 これらのデータベースに登録されているそれぞれの塩基配列には、国際的に統一されたアクセッション番号(accession number)が与えられており、 論文などでデータを引用する場合には、アクセッション番号が記載されている。 アクセッション番号は、アルファベット1文字+5桁の数字(例えば、A12345)、または、アルファベット2桁+6桁の数字(AB123456)の形式になっている。 また、塩基配列の中でタンパクをコードしている箇所(遺伝子,gene)には、タンパクの名称等に関係したニックネームが付けられている。
アクセッション番号を使って、NCBIのサイトで検索する場合は、以下のような手順で操作すればよい(データの量が大きいので、授業中に一斉に操作しないこと):
検索結果への
ショートカット。
- NCBIの検索画面にアクセスする。
- 画面上部の検索キーワードを入れる枠の中にアクセッション番号を入れる。例えば"NC000913"と入力して、searchボタンを押す。 NC000913は、分子生物学のモデル生物として重要な大腸菌(Escherichia coli)のアクセッション番号である。
- すると、かなり長大なリストが画面に表示される(ダウンロードにはそれなりに時間がかかる)。リストの先頭のLOCUS(遺伝子座)には、通常、
アクセッション番号が表示されている。次の
4639675 bp
は、 4639675塩基分のデータであることを意味する。bpは"base pair"の短縮形。 他にも、さまざまな情報が掲載されている。 - データはすでにいくつかの断片毎に仕分けられている。画面の長大なリストの中から、
/gene="lacY"
と書かれている 場所を見つけてみよう(Firefoxブラウザの「編集」メニューから「検索」)。ここで、lacYとは、ラクトース(乳糖)の分解酵素を調整する遺伝子のひとつの名称。 - 当該の行の左側に表示されている
gene
の箇所をクリックすると、lacYの塩基配列が(リストの一番おしまいに)表示される。
サンプルデータのダウンロード
NCBIに掲載されているデータと、ここで用意したサンプルデータは、
基本的には同じ遺伝子の情報であるが、対向する側の鎖の配列同士の
関係になっている。NCBIのほうのデータをよく見ると"complement(相補体)"
と書かれている行が見つかるはずだ。
演習では、(NCBIのほうではなく)、こちらのサンプルデータを用いること。
上で表示したデータは”GenBank形式”と呼ばれるフォーマットのもので、標準的な形式のひとつである。 以下で行う演習のために、GenBank形式のデータを入手し、プログラミング用の作業ディレクトリ(フォルダ)に保存しておこう。 ここでは演習のために、上記のNC000913のlacYの遺伝子のデータをあらかじめ用意しておいたので、 こちらを右クリックして、「別名でリンク先を保存」しておく。
「支援プログラム」の利用
ここでは、塩基配列のデータをプログラムで読み取って、色々な処理をすることが目標としているが、C言語についてこれまで学んだことだけを使って、
それなりに意味のあるプログラムを書くのは(無理ではないにしても)大変である。
そこで、このプロジェクトを支援するための(出来合の)Cプログラムを用意しておいた。
ここを右クリックして「別名でリンク先を保存」すると、
ga-helper.c
というファイルがダウンロードできる。先ほどの塩基配列データと合わせ、作業用のディレクトリに
ga-helper.c NC000913.gbk
の二つのファイルがあることを、まず確認する。
予行練習1:塩基配列の表示
#include <stdio.h> #include <math.h> #include "ga-helper.c" main() { int n ; int data[5000] ; n = readGBK("NC000913.gbk",data) ; printNuc(n,data) ; }
このプログラムは、GenBank形式のファイル"NC000913.gbk"を読み込んで、その中から塩基配列 に関する箇所を自動的に切り出し、配列(data)に塩基配列をセット。 次いで、その配列データの内容を、画面に表示する。
ここで、readGBK("ファイル名", 配列名)
、およびprintNuc(塩基対の数,配列名)
の箇所が、
「支援プログラム」によって提供された機能。
「支援プログラム」を使う場合には、プログラムの先頭部分に、必ず #include "ga-helper.c"
と書いておくこと。
予行演習2:データ中から特定の塩基配列を見つける
#include <stdio.h> #include <math.h> #include "ga-helper.c" main() { int i,n ; int data[5000] ; n = readGBK("NC000913.gbk",data) ; for (i=0; i<n-2; i=i+1) { if (data[i]==1 && data[i+1]==1 && data[i+2]==1) { printf("TTT was found at %d\n",i+1) ; } } }
「支援プログラム」は、塩基のTが整数値の1に, Cが2、Aが3、Gが4に、それぞれ対応するように作られている。
readGBK("NC000913.gbk",data)
を実行すると、塩基配列のデータが読み込まれ、
配列の先頭(0番目)から順に、塩基の種類に対応した数値が
3 1 4 1 3 2 1 3 1 1 1 3 3 3 3 3 3 2 3 2 3 3 3 2...
のようにセットされる。また、変数nには、塩基対の総数がセットされる。
このプログラムでは、こうして並んだデータの中に 1 (つまりT)が連続して3つ並んでいる箇所を見つけ、その位置をプリントしている。
配列の添え字は0から始まるが、「先頭から何番目」を言うときには1から順に始まるのが自然なので、(1,1,1)のパターンが見つかった場合は、
その場所として、i+1の値をプリントするようにしている。
forループの継続条件がi<n-2
になっている点も見逃せない。
現在位置から2つ先まで(計3つずつ)調べようとすると、おしまいから3つ前で止めておかないと、「行き過ぎて」しまうからだ。
予行演習3:塩基とアミノ酸の対応表を出力
#include <stdio.h> #include <math.h> #include "ga-helper.c" main() { int i,j,k ; for (i=1 ; i<=4 ; i=i+1) { for (j=1 ; j<=4 ; j=j+1) { for (k=1 ; k<=4 ; k=k+1) { printf("[%d%d%d]:",i,j,k) ; printAminoAcid(i,j,k) ; printf("\n") ; } } } }
Tは1、Cは2、Aは3、Gは4に対応しているので、
例えば、上記のプログラムで[123]
と表示される箇所は、[TCA]の塩基配列に対応する。
塩基の記号を使って表示したい場合は、プログラムの一部を
printf("[%d%d%d]:",i,j,k) ;
↓
printf("[") ; printTriplet(i,j,k) ; printf("]:") ;
と書き換えればよい。ここで、printTriplet(i,j,k);
は、数値で表現された三つの塩基の組(triplet)を、(A,T,C,G)の記号の組として表示
する「支援」機能(関数)。
また、printAminoAcid(i,j,k)
は、塩基配列[i,j,k]に対応するアミノ酸を、記号(A,..,Y)で画面に出力する関数である。
支援プログラムの機能や使い方については、こちらのページにまとめておいたので参照のこと。