ゲノムをマッピングする(リシークエンス)
- 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の割合の目安はFastQCの
Sequence 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 Recalibration
(BQSR
)をおこなうことでempirical probabilityを算出し、クオリティスコアを補正することができます。 - 既知の変異情報はvcfファイルで入力できます(『GVFとVCFを編集する』を参照)。
- まず、BaseRecalibratorで、
recalibration table
を作成し、ApplyBQSRで、bamファイル中のクオリティスコアを補正値に修正します。 -
recalibration table
は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で生成したbamファイルのQCをしておくとよいです。
- この後、たとえばゲノム全体(あるいは一部の領域)の変異をコールするには、HaplotypeCallerを用いる流れとなります。
以上。(2022年1月5日、了。)