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

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

MENU

今後のGATK解析で使用するWGSデータのマッピング -GATK解説シリーズ-part 3

FastQC, trimmomatic, bwa, samtools

導入難易度★☆☆☆☆
使用難易度★★★☆☆

今回は何をする?
前回の記事で収集したマラリア原虫のWGSデータを使ってマッピングを行い、BAM形式のファイルを生成します。内容は以前に投稿した下記の記事と重複しますが、これから解析で使うデータを扱いたいので、改めて解析を行います。

本記事では、トリミングからマッピングまでの工程をコンパクトに記載しているので、実際にご自身のデータを使って再現される場合には、この記事のスクリプトを置き替える形で使ってもらえると便利だと思います。
次回から、今回用意したデータを使ってGATK4を使用していきます。


ちなみに、以前まで使用していたWGSデータは、系統間の類似性が高すぎて使いづらかったため今後の解析には使用しない予定です。


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

以前の記事と同様に、

① 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 の検出方法を参考させていただきました。 特に、後半のスクリプト部分は非常に分かりやすく書かれていたので、ほぼ変更せず使用しています。とても参考になるサイトなので気になったらご自分で確認してみてください。

eupatho-bioinfomatics.hatenablog.com