1. 概要: bowtie2 とは
2. Mac への bowtie2 のインストール
3. インデックスの作成
4. マッピング
5. bowtie2 の結果の見方
6. bowtie2 の tips

概要: bowtie2 とは

bowtie2 とは、シークエンスから得られるリードをリファレンスに対してマッピングするためのソフトである。主に、リードが短い次世代シークエンスの結果に対して使われる。

このページに書いてあることがわかりにくい場合は、以下のリンク先をまず読んでみて下さい。

1. ターミナルの使い方
2. シェルスクリプト
3. 次世代シークエンス

Mac OSX High Sierra, 2018 年 2 月に以下の一連の操作を実行した。一部はブログへ移動。

さらに 2019 年 11 月に別の Mac でセットアップ。bowtie2, samtools ともに homebrew を使ってインストール可能だった。OS は Mojave。

Mac への bowtie2 のインストール

最初に bowtie2 のインストールを試みたときはど素人だったので、ややこしい方法を使っていた。このページには便利な方法をまとめ、かつての記事は ブログ へ移動。

まず、homebrew のページを参考にして homebrew をインストールする。

brew install bowtie2 で bowtie2 もインストール可能。確認のため which bowtie2 を実行すると、/usr/local/bin/bowtie2 となり、bin にインストールされていることがわかる。

bowtie2 の結果を解析するためには、Samtools が便利なので、これも同時にインストールしておくと良い。brew install samtools で OK。

インデックスの作成

bowtie2 index (reference, database などと呼ぶ人もいる) は fasta ファイルから bowtie2-build で作成する。たとえば sequence.fasta と使って SEQ という名前のデータベースを作成するには、以下のようにする。

bowtie2-build -f sequence.fasta SEQ

拡張子 .bt2 または .bt2l のファイルが 6 個できれば成功。この 6 個はバイナリファイルであり、中身をそのまま見ることはできない。しかし、bowtie2-inspect インデックス名で、データベースに使われている配列を見ることができる。

bowtie2 のインデックス作成に使う fasta ファイルは一続きの配列でなくても OK である。つまり、複数の配列が含まれていてもよい。ただし、それらの配列に似た部分があったりすると、read は両方にマップされることになる。bowtie2 の設定を厳しくして、そのような read を取り除くような条件でマッピングをしている場合、結果にバイアスがかかるので注意すること。

bowtie2-build のオプションの一部を表にしておく。

マッピング

コマンド bowtie2 でマッピングする。-x でインデックスを指定、-U でマップするリードを含むファイル (input) を指定する。

bowtie2 -x SEQ -1 dataset1.fastq -2 dataset2.fastq -S output.sam

インプットファイルは原則として .fastq。-f オプションをつけると .fasta も使えるらしいが、NGS の生データは一般に fastq なので、そのまま使うことが多いだろう。 圧縮ファイル fastq.gz も、解凍せずに input にできる。

以下のオプションは一部のみなので、詳細は Bowtie2 manual を参照のこと。

• リピート配列があると、マッピングはかなり影響されてしまう。これをどう取り除くかも大事。
• "*** has more quality values than read characters." というエラー。fastq ファイルは、塩基 1 個 1 個に quality value というものがついていて、そのまま解釈すると塩基数と quality value の数が一致していないためにマッピングできないというエラーである。まずはデータを再ダウンロードするのが解決への早道。大体はダウンロード中に何かが起こって、ファイルが truncated になっている。

bowtie2 の結果の見方

マッピングの結果の概要は、ターミナル上に以下のような形で表示される。数字は適当。

アウトプットのフォーマットとして、テキスト形式の .sam およびバイナリ形式の .bam がある。結果を目で見て確認したいなら、.bam よりも .sam が良い。.sam file の見方は このページ or ここ など。デフォルトのアウトプットは sam になっている気がする。

しかし、結果を引き続き解析する場合には、.bam の方がメリットが多い。一回、参考までに sam の中身を見てみて、その後は bam を基本とにするのが良いのではないかと思う。

bam でアウトプットするには、bowtie2 の結果を samtools に以下のようにパイプするのが良い。samtools のページ を参照のこと。

bowtie2 -x SEQ -U dataset.fastq | samtools view -h -bS > output.bam

bowtie2 から出力される sam file には、AS:i:5 や XS:i:5 のような項目がある。これらはalignment score であるが、どう計算されているのかマニュアルにはっきり書かれていないらしい (参考)。

Sam file には MAPQ という列があり、これはきちんと定義された値のようである。read のフィルタリングには MAPQ score を使うのがよい。ネットでざっと調べたところ、その解釈は以下の通り。不正確な記述も含まれているかもしれない。

MAPQ score は 0 - 255 の値をとる。negative, log-scaled probablity that a read is misaligned と定義される。0 は、そのリードが複数の場所にマップされたことを意味する。255 は、mapping quality not available.

2 の場合は、そのマッピングが正しくない可能性が 63%。つまり、3 つの場所にマップされうるということらしい。

bowtie2 の結果で 12924 reads (15.2%) aligned exactly 1 time という記述があり、そこから uniquely mapped read という概念が出てくる。なんとなく、複数マップされるリードはダメで、unique なものが良いという印象があり、実際にその通りなのだが、どうもこの unique という概念はあまり良くないようだ。全ては確率であり、「そのリードが正しくマップされた確率」に threshold をつけるのが妥当な計算方法。あるリードが unique に 1 箇所のみマップされたという結果が出ても、それはこの threshold 次第、つまり unique である確率も 95% とか 99% となになるので、デジタルに unique or not という形で考えない方が良い。

「bowtie2 の結果からユニークなリードを抽出する」ではなく、「MAPQ でフィルタリングして、ユニークである可能性が高いリードを抽出する」が妥当。

samtools -q 2 のようにして、MAPQ によるフィルタリングをすることができる。この場合は 2 以上のリードのみを抽出する。

このページ が良さそう。true multireads を除外したいなら、MAPQ >= 2 でよい。これは、high quality base に 4 以上のミスマッチがあるユニークリードも除外する。MAPQ 3 以上は 3 つまでミスマッチを許容。23 以上で 2 つ、40 以上で 1 つ、42 以上で 0 個。

bowtie2 の tips

Alternative splicing を調べるため、exon と intron を含むある遺伝子のゲノム配列に、transcripts をマップしたことがある。そのときの経験をメモしておく。

• 更新予定

References

1. Bioinformatics - bowtie2 (日本語) Link: Last access 2020/12/07.

コメント欄

サーバー移転のため、コメント欄は一時閉鎖中です。サイドバーから「管理人への質問」へどうぞ。

↑