ここではダウンロードしたファイルから遺伝子の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.これでアライメント作業は終了である.




TOPへ戻る