HISAT2 Samtools workflow

HISAT2→Samtoolsなワークフロー

以前書いたブログエントリなどを現状に合わせて見直して再掲載シリーズ。

HISAT2でreference genomeにmappingして、genomeに対するアラインメントを得る場合。 HISAT2のウェブサイトにすでにindexずみのそれがある場合はしなくていいが、まずはindex作成。 hisat2-buildコマンドにて。 reference genome sequenceがhogenome.fa、作成するindexの名前をhogeとすると、

# HISAT2を使うためのindex作成
% hisat2-build -p 12 hogenome.fa hoge

これはコア数が12あるMacProで動かした例(以下全て同じ)で、そこは環境に合わせて。 そして、実際のmapping。

# 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.gzfuga_2.fastq.gz、出力するSAMファイルがfuga.samとした例。 HISAT2はgzip圧縮したままのファイルでも受け付けるのと、出力はSAM形式であることに注意。

# 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が必要となるので、以下も併せて実行しておくと良い。

# BAMのindex作成
% samtools index -@ 12 fuga.bam

同じディレクトリにfuga.bam.baiというインデックスのファイルができるはず。

ちなみにBono本だと、p108-111あたりにSamoolsに関する解説がある。


Written by Hidemasa Bono in misc on 土 02 2月 2019.