Tips : [Bash][Blat] bash上で塩基配列からゲノム上の位置を特定する

　ゲノム位置の特定をネットを介さずに一瞬で済ませる方法について説明します。

Blat の導入

　まずはじめにゲノム検索に用いるBlat というプログラムを用意します。

　ダウンロードはhomebrewをインストール済みであれば、以下の命令をbashに叩いてください。
　（homebrewのインストールはこちら）

brew tap homebrew/science
brew install blat

　blatのインストールが完了したら以下のコマンドでヘルプが見れます。

blat

　この中のusageをみてください。

usage:
   blat database query [-ooc=11.ooc] output.psl

　blat ではウェブにつながず、ローカル上のファイルをデータベースとして利用します。

blatの入力

blat データベースファイル 検索する配列が入ったファイル 結果出力先ファイル

hg19.fa の取得

　では次にデータベースを用意しましょう。

　今回はヒトゲノムのリファレンスとして一般的なhg19 をデータベースとして使います。
　以下のアドレスをウェブブラウザでアクセスしてください。

hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz

　唐突に905MBのファイルがインストールされます。しばらく待ちましょう。

　ダウンロードできたら、作業をしたいディレクトリへファイルを移動します。

　以下のコマンドでファイルを展開して、hg19.faを作成します。

tar -zxvf chromFa.tar.gz;
cat chr*.fa > ./hg19.fa;

　6GBくらい食われるので注意してください。

　hg19.fa ができたら残りは捨てて構いません。

　これでデータベースファイルが取得できました。

fastaファイルにはいった配列のゲノム上の位置を配列を検索する

　今回検索する配列はこちらです。

test.fa

>1
CAGCCAACAGTGGATATTCC
>2
AGAACTACGTGGAAGTGACC
>3
GAGCTGCGCGCGGGGCCACA

ではhg19.fa があるディレクトリ上で検索しましょう。

test.faにある配列の検索

blat hg19.fa test.fa -minMatch=0 -minScore=20 output.psl

　-minScore=20 は検索配列が20塩基ということに由来します。

　これに関しては各自書き換えてください。

　しばらく待つと以下の表示がでます。

標準出力

Searched 60 bases in 3 sequences

　テキストエディタでoutput.pslは以下のようになっているはずです。

output.psl

psLayout version 3

match   mis-    rep.    N's Q gap   Q gap   T gap   T gap   strand  Q           Q       Q       Q   T           T       T       T   block   blockSizes  qStarts  tStarts
        match   match       count   bases   count   bases           name        size    start   end name        size    start   end count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
20  0   0   0   0   0   0   0   -   1   20  0   20  chr16   90354753    90128334    90128354    1   20, 0,  90128334,
20  0   0   0   0   0   0   0   +   2   20  0   20  chr1    249250621   155182200   155182220   1   20, 0,  155182200,
20  0   0   0   0   0   0   0   -   3   20  0   20  chr17   81195210    1552707 1552727 1   20, 0,  1552707,

　見辛いので拡張子を.tsvに変えてエクセルで開きます。

　T欄の
　nameが染色体番号
　start,endはゲノム上の位置となります。

　これで位置が定まりました。

　念のためtest.faの番号１の配列（CAGCCAACAGTGGATATTCC）をUCSCのgenome browserでも確認してみましょう。

　strandがマイナス（-）なので逆相補鎖が表示されていますが、間違いありませんね。

　この方法で1000個程度の配列も一気に検索できます。

　驚くべきことにその所要時間は先ほどの場合と大差ありません。

　これならblast の前で腕を組む生活にはさよならできそうです。

　また、補足として

test.faにある配列の検索２

blat hg19.fa test.fa -minMatch=0 -minScore=20 -noHead output.psl

　とすればヘッダーがない状態で保存されます。

　Rに直接データフレームとして読み込むときはこの形式がおすすめです。

　以上です。ありがとうございました。