今後のGATK解析で使用するWGSデータのマッピング -GATK解説シリーズ-part 3
FastQC, trimmomatic, bwa, samtools
導入難易度★☆☆☆☆
使用難易度★★★☆☆
今回は何をする?
前回の記事で収集したマラリア原虫のWGSデータを使ってマッピングを行い、BAM形式のファイルを生成します。内容は以前に投稿した下記の記事と重複しますが、これから解析で使うデータを扱いたいので、改めて解析を行います。
本記事では、トリミングからマッピングまでの工程をコンパクトに記載しているので、実際にご自身のデータを使って再現される場合には、この記事のスクリプトを置き替える形で使ってもらえると便利だと思います。
次回から、今回用意したデータを使ってGATK4を使用していきます。
ちなみに、以前まで使用していたWGSデータは、系統間の類似性が高すぎて使いづらかったため今後の解析には使用しない予定です。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
以前の記事と同様に、
① FastQCでデータのクオリティコントロールを行い、トリミングの閾値を設定する。
eupatho-bioinfomatics.hatenablog.com
② Trimmomaticで低クオリティーのリードを取り除くトリミングを実行する。
eupatho-bioinfomatics.hatenablog.com
③ BWA-memでリファレンスゲノムに対して、マッピングを行う。
eupatho-bioinfomatics.hatenablog.com
の3つを順番に行います。
1. 必要なツールをインストールする。
ツールのインストーツが必要な場合は、下記コマンドをターミナルに入力してインストールを実行する。
conda install -c bioconda fastqc conda install -c bioconda trimmomatic conda install -c bioconda bwa conda install -c bioconda samtools
2. 前提となるファイルについての説明。
前回入手したfastqファイルと、本ブログでいつも使用しているPlasmodium falciparum 3D7のリファレンスゲノムが必要となるので、GWASと名付けたフォルダの中にfastqと3D7_genomesというフォルダを用意して、fastqファイルとリファレンスゲノムをそれぞれ格納した。
2019年のScienceの論文で、アフリカのマラリア原虫の集団構造を2263株のWGSデータを使用して解析した大規模な研究のデータを一部拝借した。 science.sciencemag.org
フォルダの中身は以下の通りである。 (リファレンスゲノムの入手方法は別の記事で紹介しています。)
$ ls GWAS/3D7_genomes PlasmoDB-52_Pfalciparum3D7_Genome.fasta $ ls GWAS/fastq ERR015377_1.fastq.gz ERR450079_2.fastq.gz ERR015377_2.fastq.gz ERR484676_1.fastq.gz ERR015425_1.fastq.gz ERR484676_2.fastq.gz ERR015425_2.fastq.gz ERR562868_1.fastq.gz ERR063600_1.fastq.gz ERR562868_2.fastq.gz ERR063600_2.fastq.gz ERR580552_1.fastq.gz ERR1035536_1.fastq.gz ERR580552_2.fastq.gz ERR1035536_2.fastq.gz ERR636426_1.fastq.gz ERR1045266_1.fastq.gz ERR636426_2.fastq.gz ERR1045266_2.fastq.gz ERR636430_1.fastq.gz ERR1106528_1.fastq.gz ERR636430_2.fastq.gz ERR1106528_2.fastq.gz ERR701750_1.fastq.gz ERR1106529_1.fastq.gz ERR701750_2.fastq.gz ERR1106529_2.fastq.gz ERR701756_1.fastq.gz ERR2000569_1.fastq.gz ERR701756_2.fastq.gz ERR2000569_2.fastq.gz ERS032647_1.fastq.gz ERR211448_1.fastq.gz ERS032647_2.fastq.gz ERR211448_2.fastq.gz ERS032649_1.fastq.gz ERR343116_1.fastq.gz ERS032649_2.fastq.gz ERR343116_2.fastq.gz ERS347567_1.fastq.gz ERR405240_1.fastq.gz ERS347567_2.fastq.gz ERR405240_2.fastq.gz ERS347575_1.fastq.gz ERR405245_1.fastq.gz ERS347575_2.fastq.gz ERR405245_2.fastq.gz SRS399547_1.fastq.gz ERR450058_1.fastq.gz SRS399547_2.fastq.gz ERR450058_2.fastq.gz ERR450079_1.fastq.gz
3. BWAとGATK用のindexを作成する。
この後の処理に必要なインデックスファイルを作成する。
(解説シリーズ1の方法で、GATK4を導入している前提です)
bwa index -a bwtsw /3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fast samtools faidx /3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fast GATK CreateSequenceDictionary R=/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta O=/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome
上記が完了すると、Genome.fastaから下記のファイルが生成される。
ls 3D7_genomes PlasmoDB-52_Pfalciparum3D7_Genome.dict PlasmoDB-52_Pfalciparum3D7_Genome.fasta PlasmoDB-52_Pfalciparum3D7_Genome.fasta.amb PlasmoDB-52_Pfalciparum3D7_Genome.fasta.ann PlasmoDB-52_Pfalciparum3D7_Genome.fasta.bwt PlasmoDB-52_Pfalciparum3D7_Genome.fasta.fai PlasmoDB-52_Pfalciparum3D7_Genome.fasta.pac PlasmoDB-52_Pfalciparum3D7_Genome.fasta.sa
4. FastQC, trimmomatic, BWAをまとめて実行する。
各ツールの詳細は上の記事を確認してもらいたいが、今回は一度に全ての工程を実行する。
実際にはbash形式のスクリプトを作って、まとめて動かしている。
(パソコンの性能によっては非常に長時間かかることがあるので、個別に実行してください)
ステップ1 FastQCの実行。
FastQCでデータのクオリティコントロールを行い、トリミングの閾値を設定する。
大量のファイルが生成されるので、適宜フォルダに振り分ける。
workingdir=path/to/GWAS cd ${workingdir}/fastq fastqc *.fastq.gz
ステップ2 トリミング
Trimmomaticで低クオリティーのリードを取り除くトリミングを実行する。
cd ${workingdir}/fastq for fpath in `ls *_1.fastq.gz` do fname=${fpath%_1.fastq.gz} java -jar $EBROOTTRIMMOMATIC/trimmomatic-0.36.jar PE -threads 16 \ ${fname}_1.fastq.gz ${fname}_2.fastq.gz \ ${fname}_cleaned_1.fastq.gz ${fname}_unpaired_1.fastq.gz \ ${fname}_cleaned_2.fastq.gz ${fname}_unpaired_2.fastq.gz \ LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:30 done #ファイルを別のフォルダに移動する。 mkdir ${workingdir}/cleanded_fastq mv *_cleaned_*.gz ${workingdir}/cleaned_fastq/
ステップ3 マッピングを実行する。
cd ${workingdir}/cleaned_fastq for fpath in `ls *_cleaned_1.fastq.gz` do fname=${fpath%_cleaned_1.fastq.gz} bamRG="@RG\tID:L\tSM:"${fname}"\tPL:illumina\tLB:lib1\tPU:unit1" bwa mem -t 16 -R ${bamRG} $index \ ${fname}_cleaned_1.fastq.gz ${fname}_cleaned_2.fastq.gz > ${fname}.sam done #ファイルを別のフォルダに移動する。 mkdir ${workingdir}/bam mv *.sam ${workingdir}/bam/ #SAM 形式から BAM 形式に変換する。 cd ${workingdir}/bam for fpath in `ls *.sam` do samtools sort -@ 16 ${fpath} > ${fpath%.sam}.bam done
最終的なbamフォルダの中身は下のようになる。
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
お疲れ様でした。今回はこれで終わりです。
Paet3まで蛇足が長くなりましたが、次回からGATKの解説に移ります。
よければ他の記事のも見ていってください。
なお、本記事の執筆にあたりBWA および GATK4 を利用した SNPs/indels の検出方法を参考させていただきました。
特に、後半のスクリプト部分は非常に分かりやすく書かれていたので、ほぼ変更せず使用しています。とても参考になるサイトなので気になったらご自分で確認してみてください。