バイオインフォマティクスでゲノムワイド関連解析（GWAS)

バイオインフォマティクスを頑張っている方が、本ブログの内容を真似することで、自分のデータで解析ができる情報を提供することが目標です! 今はGATKの解説をメインテーマにしています。

MENU

GATK JointGenotyping -GATK解説シリーズ-part 6

GATK その他

GATK GenomicsDBimport, GATK GenotypeGVCFs, Picard VcfToIntervalList

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

eupatho-bioinfomatics.hatenablog.com

今回は何をする？

  • GATK GenomicsDBimport および GATK GenotypeGVCFs を使って、前回の記事で得たVCF形式ファイルから、変異情報を記述したローカルなデータベースを構築し、Joint Genotypingを実施して複数のvcfファイルをまとめたmerged.vcfファイルを出力します。


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

公式HPリンク
GATK GenomicsDBimport
GATK GenotypeGVCFs
Picard VcfToIntervalList

続きを読む

GATK HaplotypeCallerの使い方 後編 -GATK解説シリーズ-part 5

GATK その他

GATK HaplotypeCaller

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

eupatho-bioinfomatics.hatenablog.com

前編の記事は↓こちら eupatho-bioinfomatics.hatenablog.com

今回は何をする？

  • GATK HaplotypeCallerを使って、Part 4の記事で得たBAM形式ファイルから、変異情報の記載されたVCFファイルを出力します。
  • ものすごく時間のかかる（もし16GBのMacで実行すれば数日かかると思います）工程だったのですが、今回公式ページを確認したら、HaplotypeCallerにもマルチコアで解析するSparkモードがBeta版で搭載されていることに初めて気が付きました。開発中なのでまだ実際の解析に使用しないように注意喚起されていますが、今度試しに使ってみて感想をここに追記します。


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

公式HPのHaplotypeCallerのリンク
https://gatk.broadinstitute.org/hc/en-us/articles/360056969012-HaplotypeCaller

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

Part-4で得たBAM形式ファイルを使う。
本ブログでいつも使用しているPlasmodium falciparum 3D7のリファレンスゲノムが必要となるので、GWASと名付けたフォルダの中にbamと3D7_genomesというフォルダを用意して、markdup.bamファイルとリファレンスゲノムをそれぞれ格納した状態で開始する。（リファレンスゲノムの入手方法は別の記事で紹介しています。）

また、解説シリーズ１の方法でGATK4を導入していることを前提とする。

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

3.　Haplotype Callerを実行する。

workingdir=path/to/GWAS

cd ${workingdir}/bam
index=/path to GWAS/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta

for fpath in `ls *.markdup.bam`
do
 fname=${fpath%.markdup.bam}
 gatk HaplotypeCaller \
     --reference ${workingdir}/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta \
     --emit-ref-confidence GVCF \
     --input ${fname}.markdup.bam  \
     --output ${fname}.g.vcf \ 
     --sample-ploidy 1
done

オプションの説明
--input / -I : インプットファイル名
--output / -O: 出力ファイル名
--emit-ref-confidence / -ERC: 参照信頼度スコアを出力するモードを選択する。
--sample-ploidy / -ploidy :サンプルごとのPloidy（染色体の数）(デフォルト=2)。

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

#作成したファイルをbqsrフォルダに移動する
mkdir bqsr

mv *.g.vcf ${workingdir}/bqsr/
mv *.g.vcf.idx ${workingdir}/bqsr/


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

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

続きは↓こちら

eupatho-bioinfomatics.hatenablog.com

バイオインフォマティクス関連の書籍紹介記事 eupatho-bioinfomatics.hatenablog.com

ブログ主の自己紹介

eupatho-bioinfomatics.hatenablog.com

GATK HaplotypeCallerの使い方 前編 -GATK解説シリーズ-part 5

GATK その他

GATK HaplotypeCaller

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

eupatho-bioinfomatics.hatenablog.com

今回は何をする？

  • GATK HaplotypeCallerを使って、前回の記事で得たBAM形式ファイルから、変異情報の記載されたVCFファイルを出力します。

  • 前編では、GATK HaplotypeCallerについて勉強しましょう。


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

GATK HaplotypeCallerの概要

HaplotypeCallerは、WGSデータに含まれるSNPとインデルを、局所的なハプロタイプの再構築によって記述します。言いかえると、プログラムが変異を示す領域に遭遇すると既存のマッピング情報を破棄し、その領域を再構成します。さらに、HaplotypeCallerは、非二倍体生物のデータを扱うことができます。ただし、バリアントの可能性を計算するアルゴリズムは、（倍数性に比べて）極端な頻度の対立遺伝子には適していないので、体細胞性（癌）のバリアント検出には使用しないことが推奨されており、そのような場合には代わりにMutect2の使用が勧められています。

公式HPのHaplotypeCallerのリンク
https://gatk.broadinstitute.org/hc/en-us/articles/360056969012-HaplotypeCaller

HaplotypeCallerのワークフロー　

ワークフローのイメージ図

f:id:Harry-kun:20210718113828p:plain
Legacy GATK Forumより引用

1. アクティブな領域の定義


このステップでは、解析対象のサンプルがリファレンスと比較して変化しているゲノムの領域を選別する。これらの領域は「active regions」と定義される。バックグラウンドノイズのレベルを超す変動が見られない領域はスキップされる。この工程は、リファレンスと同一の領域で再構成を行わないことで、解析を高速化することを目的としている。

2. active regionsのアセンブリによるハプロタイプの決定


このステップでは、サンプルに存在する可能性のあるDNA配列を再構築する。active regionsにマッピングされたリードを使用して、ハプロタイプと呼ばれる配列を構築する。まず、可能性のあるハプロタイプのリストを作成するために、参照配列をテンプレートとしてactive regionsのアセンブリを構築する。

次に、各リードを順番に取得し、テンプレートとのマッチングを試みる。リードの一部がテンプレートと一致しない場合、ミスマッチを考慮してグラフに新しいノードを追加する。このプロセスを多くのリードで繰り返すと、多くの配列パターンを持つ複雑なアセンブリができあがる。しかし、GATKは各パターンを支持したリードの数を記録しているため、最も可能性の高い配列だけを選択することができる。

ハプロタイプが決定されると、それぞれのハプロタイプは元の参照配列に対して再調整され、潜在的なバリアントサイトが特定される。

3. リードデータからハプロタイプの尤度を計算する。


候補となるハプロタイプが揃ったところで、各ハプロタイプに対する各リードのペアワイズアライメントを行う。各リードを取り出し、PairHMMアルゴリズムを用いて、各ハプロタイプ（参照ハプロタイプを含む）に対して順番にアライメントを行う。PairHMMアルゴリズムでは、データの品質に関する情報（base quality scoreやindel quality scoreなど）を考慮し、各リードとハプロタイプのペアに対して、そのハプロタイプでそのリードが観測される可能性を表すスコアを出力する。

次に、これらのスコアを用いて、前のステップで特定された候補地において、個々の対立遺伝子の証拠がどれだけあるかを計算する。このプロセスはallelesに対するmarginalizationと呼ばれ、次のステップでサンプルに遺伝子型を割り当てるために最終的に使用される実際の数値を生成する。これにより、リードデータに対するハプロタイプの尤度が得られる。次に、これらの尤度データから各バリアントサイトの対立遺伝子の尤度を求める。

4. サンプルごとの遺伝子型の割り当て


前のステップでは、対象となる各候補のバリアントサイトについて、読み取りごとの対立遺伝子の尤度の表を作成した。最後に、これらの尤度を総合的に評価して、各サイトで最も可能性の高いサンプルの遺伝子型を決定する。これは、ベイズの定理を適用して、可能性のある各遺伝子型の尤度を計算し、最も可能性の高いものを選択することで行われる。これにより、ジェノタイプコールが生成されるとともに、バリアントコールが出力された場合には、出力VCFにアノテーションされる様々なメトリクスが計算される。

入力 :Variant call行うファイル (bam形式)

出力:フィルタリングされていない生のSNPおよびインデルを含むVCFまたはGVCFファイル。GVCFワークフローを利用する場合は、GVCFファイルが出力されるため、まずGenotypeGVCFsを実行し、その後フィルタリングを行ってから解析を行う必要がある。

続きは↓こちら eupatho-bioinfomatics.hatenablog.com

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

GATK その他

GATK MarkDuplicate

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

eupatho-bioinfomatics.hatenablog.com

今回は何をする？

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


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

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

続きを読む