前回の続きです。
Readをゲノムにマッピング (その１) (2017/12/19, 2018/11/19 追記あり) - Palmsonntagmorgen
今回ではbowtie, bowtie2, bwaのマッピングコマンドを説明します。
どのマッピングツールも、ゲノム配列をindex配列にまず変換し、そのindexに対してマッピングするという手順を踏みます。
直接ゲノム配列を検索するのではなく、indexを用いることで、高速・省メモリに計算することができます。
現在は、Burrows-Wheeler変換 (Burrows-Wheeler Transform; BWT) を用いてゲノムをindex化する手法が主流になっており、
上記３ツールは全てこの手法に基づいています。
Bowtieでマッピング
Bowtie, Bowtie2のindexはWebからダウンロードすることもできますが、ここではゲノム配列から自分で作成することにします。
human genome build hg38 を例とします。
index作成時のツールのバージョンとマッピング時のバージョンが異なる場合、indexの互換性が失われる可能性がありますので、
どのバージョンでindexを作成したかは記録しておきましょう。
index作成
$ bowtie-build --version # バージョン確認 $ bowtie-build hg38.fa hg38 # index作成 $ bowtie-build -C hg38.fa hg38-cs # index作成 (colorspace用)
colorspace dataを扱うことがなければ、３行目は必要ありません。
マッピング
$ bowtie -S hg38 sample.fastq -n2 -m1 -p4 > sample.hg38.sam
-S で出力フォーマットをSAMに指定しています。 -p4 は CPUを4つ使うことを示します。
-n オプションは許容ミスマッチ数を示し、-n はBLAST(MAQ)風のアラインメントスコアで評価しますが、-v オプションを使うと単純にミスマッチ数のみで評価するようになります。
-m1 はユニークマッチのみ出力するオプションで、-m10 とすると、ゲノムの10箇所以下までマップされるリードを出力します。
-k1 オプションを入れると、１リードにつきベストスコアである１箇所のみを出力します。-n2 -m1の代わりに -n2 -k1とすると、マルチリードあり、１リードにつき出力は１箇所のみ、となります。
--best --strata オプションを使った場合もマルチプルリードも出力するようになり、出力されるのは１リードにつきベストスコアである１箇所のみになります。
colorspace dataをマッピングする場合は、-Cオプションをつけ、colorspace用のインデックスを指定します。
$ bowtie -C -S hg38-cs sample.fastq -n2 -m1 -p4 > sample.hg38.sam
Bowtie2でマッピング
index作成
$ bowtie2-build --version # バージョン確認 $ bowtie2-build hg38.fa hg38 # index作成
Bowtieとほぼ同じですね。
マッピング
$ bowtie2 -x hg38 sample.fastq -p4 > sample.hg38.sam
indexは-x オプションで指定します。
Bowtie2はデフォルトの出力フォーマットがSAM形式のため、-Sオプションはありません。
出力ファイルにはマルチリードも含まれます。
マッピングのパラメータは色々あるのですが、それらの値をまとめた
--very-fast, --fast, --sensitive, --very-sensitive
というオプションが用意されています。デフォルトは--sensitiveです。
なおBowtie2は入力のreadにindelが入っている場合に分割して出力？するらしく、入力read数とSAMファイルのread数が異なる場合があります。
BWAでマッピング
index作成
$ bwa # バージョン確認 $ bwa index -p hg38 hg38.fa # index作成
BWAは何故か--versionオプションがありません。。
マッピング
$ bwa mem -t 4 hg38 sample.fastq > sample.hg38.sam
BWAもSAM形式で出力されます。出力ファイルにはマルチリードも含まれます。マッピングのパラメータは色々あるのですが、ちゃんと試したことはありません。。
PacBioとNanoporeのロングリード用のパラメータセッティングが用意されており、-xオプションで指定可能です。
その他
paired-end mappingについては記載していませんが、ほぼ同様の方法でマップすることができます。
基本的にどのツールを使っても問題ありませんが、RNA-seq, Hi-Cなどのツールの内部で特定のマッピングツールを要求されることもありますので、全て使えるようにしておくとより万全です。
