ここではダウンロードしたファイルから遺伝子のDNA配列を取り出し,さらに分析する生物群すべての配列を揃えて並べる変換(アライメント)を行う.アライメントを実行するプログラムとしてはClustalWを用いる.実行例が~shimo/class/enshu200209/にあり自由に閲覧できるのでそれを参考にすること.またこれ以降の作業全体に関する簡単なメモ(memo200209.txt)も参照.
演習を行うにあたっては
cd ~/class/enshu200209
を実行して自分のenshu200209をカレントディレクトリにすること.
1.まずNCBIのデータベースからダウンロードしたミトコンドリアゲノムのファイル(NC_001807.txt等)を~/class/enshu200209/dataディレクトリに移動する.下平が作業したディレクトリには25個のファイルを置いてある.
shimo@edu[enshu200209] cd
shimo@edu[~] cd ~/class/enshu200209
shimo@edu[enshu200209] ls data
NC_000845.txt NC_001569.txt NC_001643.txt NC_001808.txt NC_002658.txt
NC_000889.txt NC_001601.txt NC_001644.txt NC_001913.txt NC_003428.txt
NC_000891.txt NC_001602.txt NC_001645.txt NC_002008.txt NC_004023.txt
NC_001325.txt NC_001610.txt NC_001700.txt NC_002083.txt NC_004025.txt
NC_001567.txt NC_001640.txt NC_001807.txt NC_002503.txt NC_004029.txt
2.dataサブディレクトリに用意したファイルに含まれる遺伝子のDNA配列を取り出す.このためのコマンドはすでにenshu200209/binサブディレクトリにコピーしてあるcmd101である.実行するにはenshu200209ディレクトリで,
cmd101 [enter]
とする.このコマンドの中身は次のような簡単なシェルスクリプトである.
#!/usr/local/bin/bash
# setup data {catgb|nucgb|lstgb|namegb}.txt
cat data/NC_*.txt > catgb.txt
extnucgb catgb.txt > nucgb.txt
grep '#' nucgb.txt > lstgb.txt
extnamegb catgb.txt >> namegb.txt
echo "PLEASE EDIT namegb.txt"
最初の "cat data/NC_*.txt > catgb.txt" は,すべてのデータファイルを連結してenshu200209/catgb.txtというファイルにする.これから遺伝子のDNA配列を取り出してnucgb.txtというファイルにするのが"extnucgb catgb.txt > nucgb.txt"でありcmd101の作業の中心である.この実態はextnucgbというPerlスクリプトで~shimo/program/molextにあるのでファイルを覗いてみること.extnucgbの出力ファイルの実例は~shimo/class/enshu200209/nucgb.txtである.13個の遺伝子はエントリにCDSと書かれている.それ以外にも,tRNAやrRNA,さらにsourceとあるのはゲノム全体であり,これらも同時に書き出されている.#ではじまる行は入力したデータファイルの情報で,"grep '#' nucgb.txt > lstgb.txt"はこれをまとめてlstgb.txtに書き出す."extnamegb catgb.txt >> namegb.txt"はデータファイルから名前に関する情報をとりだしてnamegb.txtに書き出す.しかしこの出力は不完全なので,自分でemacsやvi等でnamegb.txtを編集して,各行が
アクセッション番号:学名:一般名
となるようにする.学名はおそらくデータファイルに含まれているが,一般名は含まれていないかもしれないので,その場合は自分で調べて書き込んでおく.実例は~shimo/class/enshu200209/namegb.txtである.
3.nucgb.txtから目的とする遺伝子だけを取り出してwork-*.nsqというファイルに書き出す.このためのコマンドは
cmd121 [enter]
である.このコマンドの中身は次のような簡単なシェルスクリプトである.
#!/usr/local/bin/bash
# setup work-*.nsq (nuc sequence data in FASTA format)
for i in ND1 ND2 COX1 COX2 ATP8 ATP6 COX3 ND3 ND4L ND4 ND5 CYTB ; do
echo $i
colseq -g \"$i\" nucgb.txt > work-$i.nsq
done
forループを使っているがもしループを使わなければ次のように12行ならべて書いても同じである.
colseq -g \"ND1\" nucgb.txt > work-ND1.nsq
colseq -g \"ND2\" nucgb.txt > work-ND2.nsq
....
colseq -g \"CYTB\" nucgb.txt > work-CYTB.nsq
ここでND1, ND2, ..., CYTBというのは遺伝子の名前である.ここでは12個の遺伝子しか解析に用いず,ND6は解析からはずしてある.(なぜか?) 出力された12個のファイルは例えばwork-ND1.nsqである.
colseqというのはPerlスクリプトで,入力したnucgb.txtの">"で始まる行に,例えば"ND1"という文字列があればそこから始まるDNAシーケンスを取り出して出力しているだけである.場合によってはこの方法がうまく行かない場合があるので,work-*.nsqファイルをemacs等で開いて中身を自分の目で確かめておくこと.
4.各遺伝子のDNA配列の書かれたファイルwork-*.nsqにコドン変換テーブル2を適用してアミノ酸配列work-*.psqを出力する.このためのコマンドは
cmd141 [enter]
である.出力サンプルはwork-ND1.psqである.
5.さてようやく配列のアライメントを行う.コマンドは
cmd201 [enter]
である.この計算には比較的時間がかかるので,バックグラウンドジョブとして
cmd201 & [enter]
として実行したほうが良いかもしれない.(あやまって何度も同じコマンドを実行しないこと! バックグラウンドジョブのことが良く理解できていないならばやらないこと.) このなかでClustalWを呼び出している.これによって,work-*.psqからアライメントファイルwork-*.alnが出力される.出力の実例はwork-ND1.alnである.実例を見て分かるように,遺伝子の長さは生物ごとに微妙に異なる場合がある.これは進化の過程でDNAシーケンスに挿入や欠損が起こったためと考えられる.アライメントでは位置を調整してギャップ(記号-)を挿入し,変化が最小になるようにつじつまを合わせている.実際にはここで用いている遺伝子ではあまり挿入や置換は起きていないが,とくに配列のはじめや最後ではギャップが挿入されている.
6.work-*.alnファイルのファイル形式を変換し,さらにギャップのある列はすべて削除してwork-*.ptnに出力する.コマンドは
cmd221 [enter]
である.ギャップのある列は後ほど行われる分析の際に悪影響を与える恐れがあるので,取り除いておく.出力の実例はwork-ND1.ptn.
7.12個のwork-*.ptnファイルをつなぎ合わせて長い配列とし,work-conc.ptnに出力する.コマンドは
cmd241 [enter]
である.出力の実例はwork-conc.ptn.これでアライメント作業は終了である.