BWAとGATKでゲノムをマッピングする

## ゲノムをマッピングする（リシークエンス）

- HiSeqXやNovaSeqの登場でIllumina系の全ゲノムショットガンリシークエンス解析が安価になりました。さらに最近はIlluminaに並ぶショートリードシークエンサーであるDNBSEQも登場しています。
- 参照ゲノム配列に対して、新しい個体の全ゲノムショットガンシークエンスデータをマッピングし、バリアント解析をするのは必須の解析ですが、あまたある方法の中でも定番のひとつである、BWA（Burrows-Wheeler Alignment Tool）からのGATK（Genome Analysis Toolkit）を使ってマッピングするためのソースをまとめました。

## Linux環境にインストールしておくべきソフト

- [BWA](http://bio-bwa.sourceforge.net/)
- [Samtools](http://www.htslib.org/)
- [GATK](https://gatk.broadinstitute.org/)（v4を使います）

## Illuminaデータの前処理

- [Illuminaデータ（fastqファイル）の前処理](https://qiita.com/hayatak/items/02e7446f54cbfc5f9259) を参照してください。
- 低品質塩基と低品質リードを除去済みの次の4つのfastqファイルを使います。
    - `hoge_fw_paired_trimmed.fastq`（forward側のペア配列）
    - `hoge_fw_unpaired_trimmed.fastq`（ペア相手のいないforward側の配列）
    - `hoge_re_paired_trimmed.fastq`（reverse側のペア配列）
    - `hoge_re_unpaired_trimmed.fastq`（ペア相手のいないreverse側の配列）

## リファレンスゲノムのindexの作成

- あらかじめ、リファレンスゲノム配列のindexを、BWA、Samtools、GATKが使えるように作成しておきます。`ref.fa`というリファレンスゲノム配列のindexを作成するのは次の書式です。

```
GENOME="ref.fa"

bwa index ${GENOME}  # BWAのためのindex化

samtools faidx ${GENOME}  # Samtoolsのためのindex化

gatk --java-options "-Xmx64G -Xms64G" \  # メモリサイズの指定
CreateSequenceDictionary -R ${GENOME}  # GATKのためのindex化
```

## BWA-MEMでのマッピング

- 定番のマッピングソフトであるBWAを使います。いくつかのアルゴリズムがありますが、ここではBWA-MEMを使います。

```
SAMPLE="hoge"
GENOME="ref.fa"

bwa mem \
-t 16 \
-R '@RG\tID:genome_hoge\tSM:hoge\tPL:ILLUMINA\tLB:hoge' \
${GENOME} \
${SAMPLE}_fw_paired_trimmed.fastq \
${SAMPLE}_re_paired_trimmed.fastq \
> ${SAMPLE}_mem_paired.sam

bwa mem \
-t 16 \
-R '@RG\tID:genome_hoge\tSM:hoge\tPL:ILLUMINA\tLB:hoge' \
${GENOME} \
${SAMPLE}_fw_unpaired_trimmed.fastq \
> ${SAMPLE}_mem_fw_unpaired.sam

bwa mem \
-t 16 \
-R '@RG\tID:genome_hoge\tSM:hoge\tPL:ILLUMINA\tLB:hoge' \
${GENOME} \
${SAMPLE}_re_unpaired_trimmed.fastq \
> ${SAMPLE}_mem_re_unpaired.sam
```

- 16スレッドを使って、それぞれペア配列、ペアのないforward側配列、ペアのないreverse側配列の3つのマッピングファイル（samファイル）を生成しました。
- `-R`で、そのリードが属するリードグループ（`@RG`）の名前（header）を定義します。これを定義しておかないと、その後のGATKでの解析ができません。`\t`で文字列を指定することができ、`ID`はその名の通りID、`SM`はサンプル名、`PL`はシークエンスのプラットフォームです。`LB`はたぶんラベルを意味しているのかと思います。（2023/11/10加筆）`-R`の中には変数が使えない（文字列として扱われる）ので、サンプル名はベタ打ちします。
- samファイルはテキスト形式ですが、その後の解析で扱いやすく、またデータサイズを小さくするために、バイナリ形式のbamファイルに変換します。そして、ソート、合体まで一気に進めます。

```
SAMPLE="hoge"

# bamファイルへ変換

samtools view -@ 16 -bS ${SAMPLE}_mem_paired.sam > ${SAMPLE}_mem_paired.bam
samtools view -@ 16 -bS ${SAMPLE}_mem_fw_unpaired.sam > ${SAMPLE}_mem_fw_unpaired.bam
samtools view -@ 16 -bS ${SAMPLE}_mem_re_unpaired.sam > ${SAMPLE}_mem_re_unpaired.bam


# ソート

samtools sort -@ 16 ${SAMPLE}_mem_paired.bam -T temp -o ${SAMPLE}_mem_paired_sorted.bam
samtools sort -@ 16 ${SAMPLE}_mem_fw_unpaired.bam -T temp -o ${SAMPLE}_mem_fw_unpaired_sorted.bam
samtools sort -@ 16 ${SAMPLE}_mem_re_unpaired.bam -T temp -o ${SAMPLE}_mem_re_unpaired_sorted.bam


# index作成

samtools index -@ 16 ${SAMPLE}_mem_paired_sorted.bam
samtools index -@ 16 ${SAMPLE}_mem_fw_unpaired_sorted.bam
samtools index -@ 16 ${SAMPLE}_mem_re_unpaired_sorted.bam


# 統計量の計算

samtools idxstats -@ 16 ${SAMPLE}_mem_paired_sorted.bam > ${SAMPLE}_mem_paired_sorted.out
samtools idxstats -@ 16 ${SAMPLE}_mem_fw_unpaired_sorted.bam > ${SAMPLE}_mem_fw_unpaired_sorted.out
samtools idxstats -@ 16 ${SAMPLE}_mem_re_unpaired_sorted.bam > ${SAMPLE}_mem_re_unpaired_sorted.out


# bamファイルの合体

samtools merge \
-@ 16 \
${SAMPLE}_mem_merged.bam \
${SAMPLE}_mem_paired_sorted.bam \
${SAMPLE}_mem_fw_unpaired_sorted.bam \
${SAMPLE}_mem_re_unpaired_sorted.bam


# ソート（合体前にソート済みなのでしなくてもいい）

samtools sort -@ 16 ${SAMPLE}_mem_merged.bam -T temp -o ${SAMPLE}_mem_merged_sorted.bam


# index作成

samtools index -@ 16 ${SAMPLE}_mem_merged_sorted.bam


# 統計量の計算

samtools idxstats -@ 16 ${SAMPLE}_mem_merged_sorted.bam > ${SAMPLE}_mem_merged_sorted.out
```

- これで`hoge_mem_merged_sorted.bam`というマッピングファイルができました。
- これはただマッピングしただけなので、いろいろ処理する必要があります。GATKで処理するのが定番です。

## PCR duplicatesの特定

- Illuminaのシークエンスは、ライブラリ作成時や、シークエンス時にPCRによる増幅をかけているため、シークエンスリードが必ずしもユニークなものではなく、いくつかのリードはPCR重複によるダブリがあります。これをPCR duplicatesといいます。PCR duplicatesの割合が多いと、それだけ偏りが生じており、正確なdepthの計算やバリアントコール（特にヘテロ接合）ができない可能性があります。そこで、PCR duplicateを特定しておき、ユニークなリードのみでの解析を行います。
- PCR duplicatesの割合の目安は[FastQC](https://qiita.com/hayatak/items/02e7446f54cbfc5f9259)の`Sequence duplication levels`として計算されているので参考になります。
- GATKの[MarkDuplicates](https://gatk.broadinstitute.org/hc/en-us/articles/360037052812-MarkDuplicates-Picard-)でduplicateのリードを特定しておきます。

```
SAMPLE="hoge"

gatk --java-options "-Xmx64G -Xms64G" \  # メモリサイズの指定
MarkDuplicates \
-I ${SAMPLE}_mem_merged_sorted.bam \
-O ${SAMPLE}_mem_merged_sorted_markdup.bam \
-M ${SAMPLE}_mem_merged_sorted_markdup_metrics.txt

# index作成と統計量計算

samtools index -@ 16 ${SAMPLE}_mem_merged_sorted_markdup.bam
samtools idxstats -@ 16 ${SAMPLE}_mem_merged_sorted_markdup.bam > ${SAMPLE}_mem_merged_sorted_markdup.out
```

## クオリティスコアの補正（BQSR）

- Illuminaシークエンスリードの塩基のクオリティスコア（base quality score）はIllumina特有の「ムラ」のようなものがあり、実態を反映しているとはいいがたい。そこで、同種ゲノムの既知の変異情報がある場合は、`Base Quality Score Recalibration`（`BQSR`）をおこなうことでempirical probabilityを算出し、クオリティスコアを補正することができます。
- 既知の変異情報はvcfファイルで入力できます（『[GVFとVCFを編集する](https://qiita.com/hayatak/items/5bb21c157a2e052d9fea)』を参照）。
- まず、[BaseRecalibrator](https://gatk.broadinstitute.org/hc/en-us/articles/360036898312-BaseRecalibrator)で、`recalibration table`を作成し、[ApplyBQSR](https://gatk.broadinstitute.org/hc/en-us/articles/360037055712-ApplyBQSR)で、bamファイル中のクオリティスコアを補正値に修正します。
- `recalibration table`は[AnalyzeCovariates](https://gatk.broadinstitute.org/hc/en-us/articles/360037066912-AnalyzeCovariates)で可視化（散布図等の作成）ができます。補正後のbamファイルをもう一度BaseRecalibratorにかけてAnalyzeCovariatesすることで、どれくらい補正されたかを可視化できます。

```
SAMPLE="hoge"
GENOME="ref.fa"
VARIANT="variant.vcf"

# BaseRecalibrator

gatk --java-options "-Xmx64G -Xms64G" \
BaseRecalibrator \
-I ${SAMPLE}_mem_merged_sorted_markdup.bam \
-R ${GENOME} \
--known-sites ${VARIANT} \
-O ${SAMPLE}_mem_merged_sorted_markdup_recal_data.table


# AnalyzeCovariates 

gatk --java-options "-Xmx64G -Xms64G" \
AnalyzeCovariates \
-bqsr ${SAMPLE}_mem_merged_sorted_markdup_recal_data.table \
-plots ${SAMPLE}_mem_merged_sorted_markdup_AnalyzeCovariates.pdf


# ApplyBQSR

gatk --java-options "-Xmx64G -Xms64G" \
ApplyBQSR \
-I ${SAMPLE}_mem_merged_sorted_markdup.bam \
-R ${GENOME} \
--bqsr-recal-file ${SAMPLE}_mem_merged_sorted_markdup_recal_data.table \
-O ${SAMPLE}_mem_merged_sorted_markdup_BQSR.bam


# BaseRecalibrator after

gatk --java-options "-Xmx64G -Xms64G" \
BaseRecalibrator \
-I ${SAMPLE}_mem_merged_sorted_markdup_BQSR.bam \
-R ${GENOME} \
--known-sites ${VARIANT} \
-O ${SAMPLE}_mem_merged_sorted_markdup_BQSR_recal_data.table


# AnalyzeCovariates after

gatk --java-options "-Xmx64G -Xms64G" \
AnalyzeCovariates \
-before ${SAMPLE}_mem_merged_sorted_markdup_recal_data.table \
-after ${SAMPLE}_mem_merged_sorted_markdup_BQSR_recal_data.table \
-plots ${SAMPLE}_mem_merged_sorted_markdup_BQSR_AnalyzeCovariates.pdf


# 統計量の計算

samtools idxstats -@ 1 ${SAMPLE}_mem_merged_sorted_markdup_BQSR.bam > ${SAMPLE}_mem_merged_sorted_markdup_BQSR.out
```

- なお、`ApplyBQSR`は自動で生成したbamファイルのindexを作成してくれます。
- これで、基本的にはマッピング（bamファイルの作成）は終了。[QualiMap](http://qualimap.conesalab.org/)で生成したbamファイルのQCをしておくとよいです。
- この後、たとえばゲノム全体（あるいは一部の領域）の変異をコールするには、[HaplotypeCaller](https://gatk.broadinstitute.org/hc/en-us/articles/360037225632-HaplotypeCaller)を用いる流れとなります。

以上。（2022年1月5日、了。）