ただ、実際には複数のソフトで結果を検証したほうが良いですし、今は当時より長い塩基数を読むことができるので、maqよりBWAの方が解析に向いているかもれない、ということで、BWAも使ってみることにしました。
まずはインストールです。
これも端末から、
sudo apt-get install BWA
で簡単にインストールできました。
あと、BWAはSNPのコールとかはやってくれないんで、そのためにsamtools
もインストールしておきます。
sudo- apt-get samtools
ですね。
では実際に使ってみたいと思います。
使うデータは前回
同様Blumenstiel et al., 2009
のデータです。
以下にまとめを書きますが、いまいちネットで情報が見つからず、試行錯誤してたどり着いた”自己流”ですので、もっと良い方法があるかもしれません。
いや、あるはずなんで、知ってたら教えてほしいです。
手順
1. ゲノムデータの索引（インデックス）を作る
bwa index 入力名(.fasta)
注: 2GB以上のデータは-a bwtswとしなくてはならない。
（defaultは-a isとなっているが、データベースの5.37倍のメモリが必要なのでここにも注意が必要）
2. fastqファイルを準備する
SRAのサイトから落として変換したfastqファイルをそのまま使うと、エラーが出てしまったので、maqのコマンドを使ってfastqファイルを分割しました。
maq fastq2bfq -n 2000000 入力名(.fastq) 出力名(.bfq)
としてファイルを分割し、bfqファイルすべてに対して以下のコマンドでfastqに戻す。
maq bfq2fastq 入力名(@XXX.bfq) 出力名(@XXX.fastq)
3. （分割した）fastqファイルをreferenceにアライメント
bwa aln 入力名(genome: @XXX.fasta) 入力名（reads: @XXX.fastq） > 出力名(@XXX.sai)
4. マッピングファイルであるSAMファイルの作成
bwa samse 入力名(genome: @XXX.fasta) 入力名（@XXX.sai） 入力名（reads: @XXX.fastq） > 出力名(@XXX.sam)
注: paired endならばsampeコマンドを使う
5. SAMファイルをBAMファイルに変換
samtools view -Sb 入力名(@XXX.sam) > 出力名(@XXX.bam)
3-5は分割したすべてのfastqファイルについて行います。
ちなみにここを自動化したので一番下にプログラムを上げておきます。
6. BAMファイルをひとつにまとめる。
samtools merge 出力名(.bam) 複数の入力ファイル名(.bam)
このBAMファイルはIGVなどのviewerで見れますし、samtoolsを使ってSNPをコールしたりできます。
SNPのコールとアノテーションについてはまた次の機会にしたいと思います。
******以下プログラムです***************************************
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
int main(void)
{
FILE
* fp;
char
row[1000], str[1000];
unsigned int i;
/*
Windows系の場合は"test.txt"にする */
fp =
fopen( "sequence.txt", "r" );
if(
fp == NULL ) {
printf( "ファイルオープンエラー\n" );
return -1;
}
while( fgets( row , sizeof( row ) , fp ) !=
NULL ) {
for(i=1; i<6; ++i){
row[strlen(row) - 1] = '\0'; /*改行と拡張子(bfq)を除く*/
}
sprintf(str,
"maq bfq2fastq bfq/%s.bfq fastq/%s.fastq", row,
row);
system(str);
sprintf(str, "bwa aln Genome/dmel-3L-chromosome-r5.32.fasta
fastq/%s.fastq > sai/%s.sai", row, row);
system(str);
sprintf(str, "bwa samse Genome/dmel-3L-chromosome-r5.32.fasta
sai/%s.sai fastq/%s.fastq > sam/%s.sam",
row, row, row);
system(str);
sprintf(str, "samtools view -Sb sam/%s.sam >
bam/%s", row, row, row);
system(str);
}
fclose( fp );
return 0;
}
**************************************************************
注: sequence.txtはbfqファイルのリストです。
