HISAT2 Samtools workflow
HISAT2→Samtoolsなワークフロー
以前書いたブログエントリなどを現状に合わせて見直して再掲載シリーズ。
HISAT2でreference genomeにmappingして、genomeに対するアラインメントを得る場合。
HISAT2のウェブサイトにすでにindexずみのそれがある場合はしなくていいが、まずはindex作成。
hisat2-buildコマンドにて。
reference genome sequenceがhogenome.fa、作成するindexの名前をhogeとすると、
#!/bin/sh
# HISAT2を使うためのindex作成
hisat2-build -p 12 hogenome.fa hoge
これはコア数が12あるMacProで動かした例（以下全て同じ）で、そこは環境に合わせて。 そして、実際のmapping。
#!/bin/sh
# HISAT2の実行
hisat2 -p 12 -x hoge -1 fuga_1.fastq.gz -2 fuga_2.fastq.gz -S fuga.sam
mapしたいpaired-end sequenceのファイルが、fuga_1.fastq.gzとfuga_2.fastq.gz、出力するSAMファイルがfuga.samとした例。
HISAT2はgzip圧縮したままのファイルでも受け付けるのと、出力はSAM形式であることに注意。
#!/bin/sh
# SAM -> BAM 変換
samtools view -@ 12 -b fuga.sam \
| samtools sort -@ 12 -T /tmp/samtools_view.$$ -o fuga.bam
そして、Samtoolsを使って、SAM→BAM変換。
同時にsortも行う例。
出てくるBAMファイルはfuga.bam。
このBAMファイルをIntegrative Genomics Viewer (IGV)などローカルなゲノムブラウザーで見る場合、indexが必要となるので、以下も併せて実行しておくと良い。
#!/bin/sh
# BAMのindex作成
samtools index -@ 12 fuga.bam
同じディレクトリにfuga.bam.baiというインデックスのファイルができるはず。
ちなみにBono本だと、p108-111あたりにSamoolsに関する解説がある。