ショートリードをゲノムにマッピングしてGATKで多型を検出する方法

GATKを素直に使うと時間が掛かりすぎるので、一部をelprepに置き換えて高速化する方法を検討したので公開しておく。
全ゲノムをいっぺんにGATKに放り込むとやはり時間が掛かるので、染色体ごとに分轄して扱う方法も併記しておく（解析のプロにとっては暗黙の了解みたいになってても、日本語ではあまり明文化されていないようだし）。

ここでは説明しないが、シーケンス解析においてSequence AlignMent (SAM) 形式のファイル、それをパソコンが読み込める形式に圧縮・翻訳したBAM（Binarized sAM）ファイル、Variant Call Format (VCF) 形式を理解していることは必須である。知らぬ者は検索せよ。

#1. 準備。まずはツールのインストール。
anacondaとmambaはここでは説明しない。必要な者は検索せよ。

conda config --add-channels bioconda --add-channels conda-forge
mamba install samtools elprep seqkit bcftools vcftools
mamba create -n gatk gatk4
mamba create -n ngsutils ngsutils

snap-aligner は https://www.microsoft.com/en-us/research/project/snap/ からダウンロードすればそのまま使える。mvで適当な場所に移してパスを通しておこう。「パスを通す」の意味を解さぬ者は検索せよ。

当たり前だが、インストールしたツールのWebサイトはいつでも参照できるようにしておくこと。特に基本的な使い方（入力ファイルと出力ファイルの記述順序など）は把握しておかないと死ぬ。

それぞれのツール用にゲノムのインデックスを作る。

snap-aligner index GENOME.FASTA GENOME.FASTA.index
seqkit seq -w 60 GENOME.FASTA > GENOME.wrapped.fasta
elprep fasta-to-elfasta GENOME.wrapped.fasta GENOME.elfasta
samtools faidx GENOME.fasta
conda activate gatk
gatk CreateSequenceDictionary -R GENOME.fasta
conda deactivate

#2. snap-alignerでショートリードをマッピング。
筆者は最近までショートリードのマッピングにはbwa-mem2を使っていたが、snap-alignerは10倍くらい高速なので乗り換えることにした。フィルタリングやbamへの変換等も全部やってくれるので圧倒的に楽になった。

snap-aligner paired GENOME.FASTA.index R1.fastq.gz R2.fastq.gz \
 -R "@RG\tID:サンプル名\tLB:library\tSM:サンプル名\tPL:Illumina"\
 -t コア数 -so -F s -F b -o -bam SAMPLE.bam

大文字の部分はあんたが実際に使うファイル名に置き換えるんやで。
-F sでMAPQ10以下のリード（複数箇所にマッピングされるリード）を除き、かつ-F bでペアの片方しかマップされなかったリードも除く。-soでファイル中のリードの並びを染色体上の位置に沿って並び替える（ポジションソートという）。

オプション：リファレンスゲノムの出来栄えにもよるが、ミスマッチを沢山含むリードが1000x以上のDepthでマッピングされる領域ができたりする。そういう領域のvariant call 処理は非常に時間が掛かる（しかも結局使えない）ので、たとえレファレンスに対してユニークにマッピングされるリードでもミスマッチが4つ以上あるものはフィルタリングしておいた方がよい。

conda activate ngsutils
samtools index SAMPLE.bam
bamutils filter SAMPLE.bam SAMPLE.filtered.bam -mismatch 3 -properpair

-mismatch はミスマッチの数を指定するフィルターだが、3つまでのミスマッチを許す、という意味なので、このオプションでミスマッチ数4以上のリードが除去される。ミスマッチのリードを除いた結果、ペアが壊れるリードも出てくるので、-properpairでフィルタリング後もペアが揃ってマップされているリードだけ残す。indexは作らなくてもbamutilsは動いてくれるが、作っておくと10倍速くらいになるので強く推奨。

#3 samtools を使って染色体ごとにbamを分割する
並列化のため。急がない人は省略可。

samtools view -bh SAMPLE.bam コンティグ名(Chr01やbctg00000000など) \
> SAMPLE.コンティグ名.bam

デキる人はループを使おう。

for i in `seq -f %02g 1 11`
do samtools view -bh SAMPLE.bam Chr$i > SAMPLE.Chr$i.bam
done

#4. elprep filter のhaplotypecallerを使ってbamからgvcfファイルを生成する
分割したBAMを使うと速い。--mark-duplicates と--remove-duplicatesでPCR duplicatesを除き、--sorting-order coordinateでファイル中のリードの並びもポジションソートをキープ。--referenceにはfastaを変換したelfastaを使う。

elprep filter SAMPLE.Chr01.bam SAMPLE.Chr01.filtered.bam \
 --mark-duplicates --remove-duplicates \
 --sorting-order coordinate \
 --reference GENOME.elfasta \
 --haplotypecaller SAMPLE.Chr01.g.vcf.gz \
 --nr-of-threads コア数

終わったら、VCFファイルのインデックスを作る。

tabix -p vcf SAMPLE.Chr01.g.vcf.gz

これらを全染色体に適用。分割したBAMでもメモリが足りない場合はelprep filterの代わりに elprep sfm ( split, filter and merge) を使う。時間は掛かるが必要なメモリは大幅に減る。

#5. GenomicsDBというデータベースを作ってそこに全サンプルのg.vcfを格納する
染色体ごとに分割しておかないとこの作業がめちゃ時間掛かる。とりあえずhaplotypecallerまでは全ゲノムでやっておいて、染色体別にGenomicsDBを作るという選択肢もあり、というか多分それが一般的（？）なやり方。

gatk --java-options “-Xmx16g” GenomicsDBImport \
 --genomicsdb-workspace-path Chr01_database \
 -R GENOME.fasta
 -V SAMPLE.Chr01.g.vcf.gz \
 -L Chr01

--java-option “-Xmx16g”はプログラムが消費するメモリ上限を上げるオプション。デフォが 4 GB だったりするので付けたほうがよいかも。上の例はメモリ上限 16 GB って意味。

複数サンプルある場合は-V SAMPLE1.Chr01.g.vcf -V SAMPLE2.Chr01.g.vdf のようにサンプルごとに-V オプションを使って指定するか、

sample1 <tab> SAMPLE1.Chr01.g.vcf.gz
sample2 <tab> SAMPLE2.Chr01.g.vcf.gz

という感じでサンプル名とファイルパスをタブ区切りで入力したテキストファイルを作って、--sample-name-map SAMPLE_LIST.txt のように指定する。

サンプル数が3桁になると1染色体分のg.vcfでも格納に数時間掛かる。というわけで、haplotypecallerを分割せずに実行した場合でも、-L で指定する染色体やコンティグを1本だけに限定して、染色体やコンティグ・スキャフォールドごとにGenomicsDBを作っておくことを強く推奨する。

なお、一度作ったGenomicsDBにサンプルを追加していくことも可能。

gatk --java-options “-Xmx16g” GenomicsDBImport \
 --genomicsdb-update-workspace-path Chr01_database \
 -V SAMPLEXX.Chr01.g.vcf.gz または –sample-name-map NewSamplesList.txt 

#6. ようやくvariant call

gatk GenotypeGVCFs -R GENOME.fasta -V gendb://Chr01_database \
-O SAMPLE1.Chr01.vcf

複数サンプルをGenomicsDBに格納したときは

gatk GenotypeGVCFs -R GENOME.fasta -V gendb://Chr01_database \
-O merged.Chr01.vcf

特定のサンプルのvcfだけが欲しいときは

gatk SelectVariants -R GENOME.fasta -V genedb://Chr01_database \
  --sample-name SAMPLE1 または -sn SAMPLE1 \
  -O SAMPLE1.Chr01.vcf

#7. 信頼性の高いSNPやINDELだけをフィルタリングする
まずはVCFからSNPとINDELを分離する。

SNPサイトの抽出

gatk SelectVariants -R GENOME.fasta -V SAMPLE.Chr01.vcf または merged.vcf \
  -O SAMPLE1.Chr01.snp.vcf または merged.snp.vcf \
  --restrict-alleles-to BIALLELIC \
  --select-type-to-include SNP

BIALLELICは対立遺伝子が2種類しかない、つまり変異型のSNPが1種類しかないものだけを抽出するの意。multi-allelicなサイトも集団遺伝学的には意味があるが、解析が複雑化するので殆どの場合無視される。

INDELサイトの抽出

gatk SelectVariants -R GENOME.fasta -V SAMPLE.Chr01.vcf または merged.vcf \
  -O SAMPLE.Chr01.indel.vcf または merged.indel.vcf \
  --select-type-to-include INDEL

それぞれから信頼性の高い多型サイトのみを抽出する。
以下の処理では大量のWarningメッセージが出力されるが、所詮警告であってエラーではない、ということで気にしてはいけないらしい（早く直して欲しいなぁ）。

SNPのフィルタリング

gatk VariantFiltration -R GENOME.fasta \
  -V SAMPLE.Chr01.snp.vcf または merged.snp.vcf \
  -O SAMPLE.Chr01.snp.filtered.vcf または merged.snp.filtered.vcf \
  -filter "QD < 2.0" --filter-name "QD2"       \
  -filter "QUAL < 30.0" --filter-name "QUAL30" \
  -filter "SOR > 4.0" --filter-name "SOR4"     \
  -filter "FS > 60.0" --filter-name "FS60"     \
  -filter "MQ < 40.0" --filter-name "MQ40"     \
  -filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
  -filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8"

-V が入力するVCFで-Oが出力するVCF。あとはフィルタリングに用いる指標とその閾値、そしてフィルタ名 MQはroot mean square Mapping Quality, QDはQuality by Depth, FS はFisher Strand、MQRankSumはMapping Quality Rank Sum Test, ReadPosRankSum は Read Position Rank Sum Test の略。詳しくは↓を参照のこと。
https://gatk.broadinstitute.org/hc/en-us/articles/360035890471-Hard-filtering-germline-short-variants

INDELのフィルタリング

gatk VariantFiltration -R GENOME.fasta \
  -V SAMPLE.Chr01.indel.vcf または merged.indel.vcf \
  -O SAMPLE.Chr01.indel.vcf または merged.indel.filtered.vcf \
  -filter "QD < 2.0" --filter-name "QD2"       \
  -filter "QUAL < 30.0" --filter-name "QUAL30" \
  -filter "FS > 200.0" --filter-name "FS200"   \
  -filter "SOR > 10.0" -filter-name "SOR10"    \
  -filter "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20"

ポジションに関する閾値が少し甘くなっています。まあINDELだからね。マッピングの位置がズレてしまうことがあるわけで（どういうことか考えてみよう）。

終わったら、フィルタリング前後でどこが変わっているか、less -S を使って中身を見比べてみること。

フィルタリングにPASSした多型サイトだけを抽出しておくと後々便利。ファイルも小さくなるし。

bcftools view -f PASS -Oz -o SAMPLE.Chr01.snp.pass.vcf.gz SAMPLE.Chr01.snp.filtered.vcf
bcftools view -f PASS -Oz -o SAMPLE.Chr01.indel.pass.vcf.gz  SAMPLE.Chr01.indel.filtered.vcf

#8. Base Quality Score Recalibration
フィルタリングを終えたVCFをもとにBase Quality Score Recalibration (BQSR) を行って、リードのQscoreを修正する。リード側のエラーに惑わされにくくするということかな。詳しくは”gatk bqsr”で検索せよ。

まずは仮のvariant callから多型サイトを「既知の多型サイトの一覧表」へと変換する。

elprep vcf-to-elsites SAMPLE.Chr01.snp.pass.vcf.gz SAMPLE.Chr01.snp.elsites
elprep vcf-to-elsites SAMPLE.Chr01.indel.pass.vcf.gz SAMPLE.Chr01.indel.elsites

merged.Chr01.vcfの場合はそれを使う。

elprep vcf-to-elsites merged.Chr01.snp.pass.vcf.gz merged.Chr01.snp.elsites 
elprep vcf-to-elsites merged.Chr01.indel.pass.vcf.gz merged.Chr01.indel.elsites

そして再びelprep filter のhaplotypecaller

elprep filter SAMPLE.Chr01.filtered.bam SAMPLE.Chr01.filtered.2.bam \
 --mark-duplicates --remove-duplicates \
 --sorting-order coordinate \
 --reference GENOME.elfasta \
 --bqsr SAMPLE1.recal
 --known-sites SAMPLE.Chr01.snp.elsites または merged.Chr01.snp.elsites \
 --known-sites SAMPLE.Chr01.indel.elsites または merged.Chr01.indel.elsites \
 --haplotypecaller SAMPLE.Chr01.bqsr.g.vcf.gz \
 --nr-of-threads コア数

tabix -p vcf SAMPLE.Chr01.bqsr.g.vcf.gz

#9. haplotypecaller ~ VariantFiltration をもう一回
BQSRが終わったら、5．GenomicsDBImport と 6．GenotypeGVCFsと 7．SelectVariants & VariantFiltration を全く同じ手順で繰り返すと、「順当な手続きを踏んだ」 Variant Call Format（VCF）ファイルが得られる。お疲れ様。