BWA(ヒト single-end)

BWAを利用したマッピングの流れ

BWA は短いリードのみならず、長いリードのマッピングにも対応している。BWA を利用したマッピングは左図のように、インデックスの作成とマッピングの 2 ステップからなる。

このページでは、SRR1163657 と SRR1163658 の 2 つのサンプルデータを用いて、BWA の使い方を示す。リファレンスゲノムは Ensembl からヒトゲノムをダウンロードして利用する。(※サンプルデータはクオリティコントロール後のデータとする)

  1. リファレンスゲノムの取得
  2. インデックスの作成
  3. マッピング
  4. カウントデータの取得
 

リファレンスゲノムの取得

ヒトの全ゲノムデータは公共データベースである Ensembl でダウンローできる。しかし、全ゲノムをダウンロードしても、どの領域が遺伝子でどの領域がそうではないかのアノテーション情報が含まれていない。そのため、(マッピングには必要はないが、)全ゲノムの他にアノテーションデータも合わせてダウンロードすると便利である。Ensembl のダウンロードページにアクセスすると以下の様な項目が見られる。

BWAによるマッピング

Ensembl のダウンロードページにて、上図の赤い枠で囲まれたファイルをダウンロードする。

該当リンクをクリックした後に、染色体別のデータなど様々なタイプのデータが表示される。ここでダウンロードするのは全ゲノムデータであるので、以下のファイル名のファイルをダウンロードする。

  • Homo_sapiens.GRCh38.dna.toplevel.fa.gz

また、同時にアノテーションファイル(GTF)もダウンロードしておく。

  • Homo_sapiens.GRCh38.gtf.gz

ダウンロードしたファイルは gzip で圧縮されているため、ターミナルから以下のコマンドを実行して解凍させる。

gunzip Homo_sapiens.GRCh38.dna.toplevel.fa.gz
gunzip Homo_sapiens.GRCh38.gtf.gz

展開後、ディレクトリの中には Homo_sapiens.GRCh38.dna.toplevel.fa と Homo_sapiens.GRCh38.gtf の 2 つのファイルができる。

インデックスの作成

マッピングに先立ち、ダウンロードしたリファレンスゲノムに対してインデックスを作成する。インデックスの作成は bwa index のコマンドを実行して作成する。オプション -p はインデックスに名前を付けるというオプションである。

以下の例では、Homo_sapiens.GRCh38.dna.toplevel.fa のインデックスを作成し、そのインデックスの名前を human_index としている。

bwa index -p human_index Homo_sapiens.GRCh38.dna.toplevel.fa

ヒトゲノムの場合、作成時間が、パソコンの性能にもよるが、 20 時間前後を要する。

マッピング

BWA には 3 つのアルゴリズムが実装されている。そのなかで、著者らは BWA-MEM を推奨している。ここでは BWA-MEM を利用してマッピングを行う。

ここで利用するサンプルデータは SRR1163657 と SRR1163658 である。両サンプルデータともに、5' 末端の ACGT の割合が不自然で、3' 末端のクオリティが低いため、両端をトリムした。(参照:クオリティコントロール

マッピングは以下のように行う。

bwa mem human_index SRR1163657.fastq > SRR1163657.sam
bwa mem human_index SRR1163658.fastq > SRR1163658.sam

マッピング結果は SAM 形式のファイルに保存される。複数の CPU でマッピングを行う場合は -t オプションを利用する。

bwa mem -t 8 human_index SRR1163657.fastq > SRR1163657.sam
bwa mem -t 8 human_index SRR1163658.fastq > SRR1163658.sam

カウントデータの取得

カウントデータの取得は HTSeq 中の htseq-count コマンドを利用する。この際に、マッピング結果の SAM ファイルとともに、アノテーションファイル(GTF)を利用する。

htseq-count SRR1163657.sam Homo_sapiens.GRCh38.gtf > SRR1163657.txt
htseq-count SRR1163658.sam Homo_sapiens.GRCh38.gtf > SRR1163658.txt

カウントデータはそれぞれ SRR1163657.txt と SRR1163658.txt に保存される。

カウントデータが複数のファイルにまたがっているため、これらを join コマンドで結合する。また、上述のカウントデータが記載されいてるファイルのなかに、遺伝子のカウントデータの他に、アノテーションがないリード数やマップされてなかったリード数なども保存されている。前者の場合は ENSG から始まり、後者の場合は「__no_feature」や「__not_aligned」などとなっている。そのため、grep 文を利用して ENSG だけから始まるデータのみを取得する。

# カウントデータのみを tmp.count.txt (一時ファイル)に保存する
join -j 1 SRR1163657.count.txt SRR11693658.count.txt > tmp.count.txt

# ヘッダーを all.count.txt に書き入れる
echo "EnsemblID    SRR1163657    SRR1163658" > all.count.txt

# カウントデータを all.count.txt に書き入れる
grep "^ENSG" tmp.count.txt >> all.count.txt

# 一時ファイルを削除する
rm tmp.count.txt

References

  1. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009, 25(14):1754-60. PubMed Abstract
  2. Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010, 26(5):589-95. PubMed Abstract