Maqの使い方まとめ

自分用まとめです

下調べ

GCOEテクニカルセミナーのレジメがすごく役立った(ここ)
S.pombe分裂中さんの記事もわかりやすかったです(ここ)
参考にさせていただきました ありがとうございます

レファレンスのシークエンスをとってくる

Humanのレファレンス配列をUCSCからとることにする
UCSC Genome Browser: DownloadsのHuman Genome (hg19, GRCh37)
Data set by chromosomeにある
染色体毎の配列(chr1.fa.gzなど)をダウンロード


解凍して、一つのファイルにまとめておく

cat chr*.fa > allchr.fa

シークエンサーのデータをとってくる

ウチには次世代シークエンサーは無いので(^^;)、NCBISRA(Short Read Archive)からとってくる

今回の目的に合いそうなものとして
Study: Targeted capture and massively parallel sequencing of human exomes JS0001 (SRP000910)
のデータを使うことにする


そのうちの一つ
Human HapMap individual NA19240
のデータ(SRX005930)のSRR017992をダウンロードしてみた。解凍してみたら10GB近くあった(@_@)


ところで、
SRAのFastqデータをMaqで使うには、ちょっと加工してやらなければならないらしい
タイトル行にあるスペースを削除すればいいそうな GCOEのテキストに習ってsedで加工してやる

sed 's/ //g' < SRR017992.fastq > SRR017992.fastq_ns

MAQでマッピング

いろいろ細かいことは後でマニュアルを読んで調べることとして、まずはeary runでマッピングしてみた

$ maq.pl earyruy -1 60 -d SRR017992 allchr.fa SRR017992.fasta_ns

-1 60 のオプションは、Fastq配列の頭から60ベースのみをマッピングに使うためにつけた*1


後は待つこと半日〜1日
ウチのマシン(2x2.26 GHz Quad Core Xeon, 8GB)で36,394,026spotsのデータをマッピングするのに、約24時間かかった

マルチコアに最適化していないと作者もいっているとおりで、シングルコアしか使っていなかった。


MaqViewrでデータを眺める

まずはインデックスをつけてやる必要がある

$ maqindex -i -c consensus.cns all.map


そしたらMaqviewでデータを眺めてみる

$ maqview -c consensus.cns all.map

*1:このfastqデータは72ベースあるのだが、おしりの部分のクオリティが低くマッピングから除外した。そうしないとMaqがエラーで止まってしまう