4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

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

Last updated at Posted at 2022-01-04

ゲノムをマッピングする(リシークエンス)

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

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

Illuminaデータの前処理

  • Illuminaデータ(fastqファイル)の前処理 を参照してください。
  • 低品質塩基と低品質リードを除去済みの次の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の割合の目安はFastQCSequence duplication levelsとして計算されているので参考になります。
  • GATKのMarkDuplicatesで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 RecalibrationBQSR)をおこなうことでempirical probabilityを算出し、クオリティスコアを補正することができます。
  • 既知の変異情報はvcfファイルで入力できます(『GVFとVCFを編集する』を参照)。
  • まず、BaseRecalibratorで、recalibration tableを作成し、ApplyBQSRで、bamファイル中のクオリティスコアを補正値に修正します。
  • recalibration tableAnalyzeCovariatesで可視化(散布図等の作成)ができます。補正後の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で生成したbamファイルのQCをしておくとよいです。
  • この後、たとえばゲノム全体(あるいは一部の領域)の変異をコールするには、HaplotypeCallerを用いる流れとなります。

以上。(2022年1月5日、了。)

4
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?