GATK MarkDuplicatesの使い方 -GATK解説シリーズ-part 4

GATK MarkDuplicate

使用難易度★★☆☆☆
本記事は、GATK解説シリーズのPart 4です。

eupatho-bioinfomatics.hatenablog.com

今回は何をする？

• 前回の記事で得たBAM形式ファイルを使って、GATK MarkDuplicates/MarkDuplicateSparkにより重複したリードにタグを付けます。

Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。 フォローする @harrykun_blog

1. GATK MarkDupicatesの概要

このツールでは、DNAライブラリ中の単一DNA断片に由来する重複リードを検出してタグを付けることができます。

BAMファイルまたはSAMファイル内の重複リードを検出してタグ付けします。重複リードとは、単一のDNA断片に由来するリードと定義され、PCRを使用したライブラリ構築などのサンプル調製時に発生する可能性があります。また、PCRの重複とは別に光学的重複も検出します。光学的重複とは、1つの増幅クラスターがシーケンシング装置の光学センサーによって複数のクラスターとして誤検出されることによって生じる重複です。

PCRによる重複の検出は、EstimateLibraryComplexityに記載の方法に従って行われ、BAM形式ファイル内のペアリードの配列から重複を判定します。各リードのN塩基から5塩基（デフォルト）を使ってソートされ、ペア同士が互いにギャップなく一致し、全体のミスマッチ率がMAX_DIFF_RATE（デフォルトでは0.03）以下であれば重複しているとみなされます。重複リードを収集した後、base-quality scores でリードをランク付けするアルゴリズムを使用して、プライマルリードと重複リードを区別します。このツールの出力は、各リードの重複が識別された新しいBAMファイルです（重複リードを破棄するわけではない）。

以前は、重複かどうかがだけ判別されましたが、重複の種類を特定するものではありませんでした。このため、最近BAMファイルのオプション出力として、DT（duplicate type）という新しいタグが追加されました。TAGGING_POLICYオプションを使用すると、すべての重複をマークする（All）、光学的重複のみをマークする（OpticalOnly）、または重複をマークしない（DontTag）ようにプログラムに指示することができます。BAMファイルの出力は、ライブラリPCRで生成された重複（LB）、またはシーケンシングプラットフォームのアーティファクトの重複（SQ）のいずれかのDT値を持ちます。必要に応じてREMOVE_DUPLICATEまたはREMOVE_SEQUENCING_DUPLICATESを使うことで、重複リードを削除することができます。

2. MarkDuplicatesSparkとMarkDuplicatesの違い

細かな違いはいくつかあるようですが、簡単にいうと複数コアで実行可能になったものがSparkだと認識しています（間違っているかも？）。16GB以上のメモリを持つPCでの使用が推奨されているようです。時間があったら速度を比較してみて、ここに追記したいと思います。

公式HPのリンク

MarkDuplicates
https://gatk.broadinstitute.org/hc/en-us/articles/360057439771-MarkDuplicates-Picard-

MarkDuplicatesSpark
https://gatk.broadinstitute.org/hc/en-us/articles/360057438771-MarkDuplicatesSpark

3.　前提となるファイルについての説明。

前回得たBAM形式ファイルを使う。 また、解説シリーズ１の方法でGATK4を導入していることを前提とする。

2019年のScienceの論文で、アフリカのマラリア原虫の集団構造を2263株のWGSデータを使用して解析した大規模な研究のデータを一部拝借した。 science.sciencemag.org

GWASフォルダの中に、bamフォルダを置き、以下のbamファイルを格納した状態で開始する。

#workingdir=path/to/GWAS
ls ${workingdir}/bam

ERR015377.bam        ERR484676.bam
ERR015425.bam        ERR562868.bam
ERR1035536.bam       ERR580552.bam
ERR1045266.bam       ERR636426.bam
ERR1106528.bam       ERR636430.bam
ERR1106529.bam       ERR701750.bam
ERR2000569.bam     ERR701756.bam
ERR211448.bam        ERS032647.bam
ERR343116.bam       ERS032649.bam
ERR405240.bam      ERS347567.bam
ERR405245.bam      ERS347575.bam
ERR450058.bam      SRS399547.bam
ERR450079.bam   

4.　MarkDuplicatesを実行する。

cd ${workingdir}/bam
for fpath in `ls *.bam`
do
 gatk MarkDuplicates　\
            --INPUT ${fpath%.bam}.bam \
            --OUTPUT ${fpath%.bam}.markdup.bam \
            -M ${fpath%.bam}.markdup.metrics.txt            
done

オプションの説明
--INPUT / -I : インプットファイル名（マッピングされていなくても可）
--OUTPUT / -O: 出力ファイル名
--metrics-file / -M :メトリクスファイルを出力する際のファイル名

補足 MarkDuplicatesSparkを実行する。

cd ${workingdir}/bam
for fpath in `ls *.bam`
do
 gatk MarkDuplicatesSpark \
            --INPUT ${fpath%.bam}.bam \
            --OUTPUT ${fpath%.bam}.markdup.bam \
            -M ${fpath%.bam}.markdup.metrics.txt \
            --conf 'spark.executor.cores=16'
done

オプションの説明
--conf 'spark.executor.cores=16': コア数

解析が終了するとbamフォルダにファイルが生成される。

XXXX.markdup.bam
XXXX.markdup.bam.bai
XXXX.markdup.bam.sbi
XXXX.metrics.txt

お疲れ様でした。今回はこれで終わりです。
よければ他の記事のも見ていってください。

なお、本記事の執筆にあたりBWA および GATK4 を利用した SNPs/indels の検出方法を参考させていただきました。 とても参考になるサイトなので気になったらご自分で確認してみてください。

バイオインフォマティクス関連の書籍紹介は↓こちら eupatho-bioinfomatics.hatenablog.com