GATK4をsplit intervalを使って高速化する

TL;DR

GATK4は実行に時間がかかるツールですが、マシンパワーさえあればsplit intervalを使って高速化できます。split intervalについては日本語文献が見つからなかったのでまとめておきます。Sparkを使った実装も進められているようでうが、まだ全てに対応しているわけではないようです(2020/09/28現在)。

split intervalディレクトリの作成

まずpicard ScatterIntervalsByNsを使ってinterval listを作成します。リファレンスゲノムのからインターバルリスト形式のファイルに変換します。ポジションとATGC、Nの数で構成されたファイルです。

picard ScatterIntervalsByNs REFERENCE=ref.fa OUTPUT_TYPE=ACGT OUTPUT=ref.interval_list

次にsplit intervalを作成します。インターバルリストの分割です。最終的にすべての分割されたインターバルは同一の塩基数を持ちます。gatk4 SplitIntervalsを使います。この例では12個に分割しています。

gatk4 SplitIntervals -R ref.fa -L ref.interval_list --scatter-count 12 -O interval_list_12

split intervalつきのGATK4の実行例

BQSR -> ApplyBQSR -> HaplotypeCallerくらいの実行例を載せておきます。基本的には、-Lオプションでインターバルリストを指定する、forで分割されたものについて回す、という感じです。あとwaitでちゃんと処理終了を待つ必要があります。

最終的に分割されたvcfをpicard GatherVcfsで集めています。

ref=/path/to/fasta
ref_dir=/path/to/interval_list_12
basename=your_file_basename

for i in `seq -f '%04g' 0 11`; do
    outfile=${basename}_recal_data_$i.table
    gatk --java-options "-Xmx4G" \
        BaseRecalibrator \
        -L $ref_dir/interval_list_12/${i}-scattered.interval_list \
        -R $ref \
        -I ${basename}.bam \
        --known-sites /path/to/known.vcf \
        -O $outfile &
done
wait

for i in `seq -f '%04g' 0 11`; do
    bqfile=${basename}_recal_data_$i.table
    output=${basename}_recal_$i.bam
    gatk --java-options "-Xmx4G" \
        ApplyBQSR \
        -L $ref_dir/interval_list_12/${i}-scattered.interval_list \
        -R $ref \
        -I ${basename}.bam \
        -bqsr $bqfile \
        -O $output &
done
wait

for i in `seq -f '%04g' 0 11`; do
    infile=${basename}_recal_$i.bam
    outfile=${basename}_$i.g.vcf
    gatk --java-options "-Xmx4G" \
        HaplotypeCaller \
        -L $ref_dir/interval_list_12/${i}-scattered.interval_list \
        -R $ref \
        -I $infile \
        -ERC GVCF \
        -O $outfile &
done
wait

for i in `seq -f '%04g' 0 11`; do
    infile=${basename}_$i.g.vcf
    gatk --java-options "-Xmx4G" \
        GenotypeGVCFs \
        -L $ref_dir/interval_list_12/${i}-scattered.interval_list \
        -R $ref \
        -V $infile \
        -O ${basename}_$i.vcf &
done
wait

in_vcfs=`echo $(ls ${basename}_00*.vcf | grep -v ".g.vcf") | sed -e 's/ / I=/g'`

picard \
    GatherVcfs \
    R=$ref \
    I=$in_vcfs \
    O=$basename.vcf

参考

この記事に関するIssueをGithubで作成する

Read Next