Pythonプログラミング（総合演習課題）

このページでは、塩基配列データをCのプログラムで扱うための具体的な手順を説明する。

プログラミングのための準備

塩基配列データの入手

現在までに知られている塩基配列の多くは、データベースを使って簡単に検索することができる。例えば、 Genbank(NCBI)のヌクレオチド検索データベース、日本DNAデータバンク(DDBJ)の検索ページ、 EMBL-EBIのデータベース、などは良く知られている。それらの中から、必要なデータを絞り込むには、専門的な知識や用語を知っておく必要があるので、こうしたデータベースを自在に使いこなせるようになには、専門の勉強を積む必要があるだろう（何せ、この講義の担当教員の専門は物理なので、「教養」程度のことしか分からない)。

そうは言っても、検索自体は、簡単に行うことができる。これらのデータベースに登録されているそれぞれの塩基配列には、国際的に統一されたアクセッション番号(accession number)が与えられており、論文などでデータを引用する場合には、アクセッション番号が記載されている。アクセッション番号は、アルファベット1文字+5桁の数字（例えば、A12345)、または、アルファベット２桁+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の遺伝子のデータをあらかじめ用意しておいたので、こちらを右クリックして、「別名でリンク先を保存」しておく。

「支援プログラム」の利用

支援プログラムのダウンロード

ここでは、塩基配列のデータをプログラムで読み取って、色々な処理をすることが目標としているが、Pythonについてこれまで学んだことだけを使って、それなりに意味のあるプログラムを書くのは（無理ではないにしても）大変である。そこで、このプロジェクトを支援するための（出来合の）コードを用意しておいた。ここを右クリックして「別名でリンク先を保存」すると、 ga.pyというファイルがダウンロードできる。先ほどの塩基配列データと合わせ、作業用のディレクトリに

ga.py       NC000913.gbk

の二つのファイルがあることを、まず確認する。

予行練習１：塩基配列の表示

# coding: utf-8

import ga

data = ga.readGBK("NC000913.gbk")
ga.printNuc(data)

このプログラムは、GenBank形式のファイル"NC000913.gbk"を読み込んで、その中から塩基配列に関する箇所を自動的に切り出し、配列（data）に塩基配列をセット。次いで、その配列データの内容を、画面に表示する。

ここで、data=ga.readGBK("ファイル名")、およびga.printNuc(リスト名)の箇所が、「支援プログラム」によって提供された機能。「支援プログラム」を使う場合には、プログラムの先頭部分に、必ず import gaと書いておくこと。

予行演習２：データ中から特定の塩基配列を見つける

# coding: utf-8
import ga

data = ga.readGBK("NC000913.gbk")

for i in range(len(data)-2):
    if data[i]==1 and data[i+1]==1 and data[i+2]==1:
        print("TTT was found at",i+1)

「支援プログラム」は、塩基のTが整数値の1に, Cが２、Aが3、Gが4に、それぞれ対応するように作られている。 ga.readGBK("NC000913.gbk")を実行すると、塩基配列のデータが読み込まれ、リストの先頭（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...

のようにセットされる。

このプログラムでは、こうして並んだデータの中に 1 （つまりT）が連続して３つ並んでいる箇所を見つけ、その位置をプリントしている。配列の添え字は０から始まるが、「先頭から何番目」を言うときには１から順に始まるのが自然なので、(1,1,1）のパターンが見つかった場合は、その場所として、i+1の値をプリントするようにしている。 forループの範囲がrange(len(data)-2) になっている点も見逃せない。現在位置から２つ先まで（計３つずつ）調べようとすると、おしまいから３つ前で止めておかないと、「行き過ぎて」しまうからだ。

予行演習３：塩基とアミノ酸の対応表を出力

# coding: utf-8
import ga

for i in range(1,5,1):
    for j in range(1,5,1):
        for k in range(1,5,1):
            print("[",i,j,k,"]:",sep="",end="")
            ga.printAminoAcid(i,j,k)
            print("")

Tは1、Cは２、Aは3、Gは4に対応しているので、例えば、上記のプログラムで[123]と表示される箇所は、[TCA]の塩基配列に対応する。塩基の記号を使って表示したい場合は、プログラムの一部を


print("[",i,j,k,"]:",sep="",end="")
　　　↓
print("[", ga.printTriplet(i,j,k), "]:",sep="",end="")

と書き換えればよい。ここで、ga.printTriplet(i,j,k)は、数値で表現された三つの塩基の組(triplet)を、(A,T,C,G)の記号の組として表示する「支援」機能（関数）。

また、printAminoAcid(i,j,k) は、塩基配列[i,j,k]に対応するアミノ酸を、記号（A,..,Y）で画面に出力する関数である。

支援プログラムの機能や使い方については、こちらのページにまとめておいたので参照のこと。