SNPs や indels などの short variants を検出するプログラムの中で、とりわけ GATK がよく使われている。解析の流れとして、基本的に BWA または STAR でマッピングを行い、そのマッピング結果を GATK で解析し variant calling を行う。以下に、ニューヨーク大学の Genomics Core ウェブサイトで公開されている解析パイプラインと GATK の Best Practices Workflows を参考にして作成した解析パイプラインを示す。以下の解析パイプラインを大まかな解析の流れとして捉えて、必要に応じて個々の処理をカスタマイズして使うとよい。
解析環境の準備
この解析では、QC、マッピングや variant calling など多数操作を行う必要がある。それぞれの操作で途中結果が生成されるため、これらの結果が混同しないように、予め解析ディレクトリおよびフォルダ構造をここで決めておくよい。そこで、以下のように、atgwas フォルダを用意して、その下に様々なフォルダを作成する。また、atgwas フォルダを中心に操作するため、このフォルダを変数として宣言し、あとで使いやすいようにする。
# project path
PROJECT_PATH=/home/biopapyrus/projects/atgwas
mkdir -p ${PROJECT_PATH}
cd ${PROJECT_PATH}
pwd
## /home/biopapyrus/projects/atgwas
mkdir genome
mkdir fastq
mkdir cleaned_fastq
mkdir bam
mkdir bqsr
mkdir vcfデータセット
ここでは Arabidopsis thaliana 1001 Genomes プロジェクト のデータの一部をサンプルとして使用する。1001 系統のデータのうち、今回は SRR094979、SRR095680、SRR095786、SRR095851、SRR095698 系統のデータのみを使用する。これらのデータは ENA からダウンロードできる。なお、オリジナルのデータサイズが大きく、処理時間が長い。テスト実行であれば、この時点で FASTQ のサイズを減らすと良い。
cd ${PROJECT_PATH}/fastq
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR094/SRR094979/SRR094979_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR094/SRR094979/SRR094979_2.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095680/SRR095680_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095680/SRR095680_2.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095786/SRR095786_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095786/SRR095786_2.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095851/SRR095851_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095851/SRR095851_2.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095698/SRR095698_1.fastq.gz
wget ftp.sra.ebi.ac.uk/vol1/fastq/SRR095/SRR095698/SRR095698_2.fastq.gz次に、リファレンスとなる配列を Ensembl からダウンロードする。ここでは Arabidopsis thaliana の最新のゲノム配列 TAIR10 を使用する。ゲノム配列とアノテーションを同時にダウンロードする。ダウンロード後にこれらのファイルを展開する。
cd ${PROJECT_PATH}/genome
# reference sequence
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# gene annotation
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-40/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.40.gff3.gz
# decompress
gunzip Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.40.gff3.gzQC
FASTQ ファイルの中のリードに対して、クオリティの低いものをトリムしたり、削除したりするクオリティコントロールを行う。まず、ダウンロードした直後のデータのクオリティを確認するために、QC レポートを作成する。
cd ${PROJECT_PATH}/fastq
fastqc *.gzQC レポートを参考にしてクオリティコントロールの閾値を設定する。その後、Trimmomatic でフィルタリングを行う。なお、ここで使用したサンプルデータは 10 年前のデータであるために、クオリティが非常に悪い。厳しくフィルタリングするとリードがなくなるので、ここでは非常に緩い閾値を設定してフィルタリングを行なっている。最近のシーケンサーから出力されたデータを使用する場合は、閾値を厳しく設定しおくと良い。
cd ${PROJECT_PATH}/fastq
for fpath in `ls *_1.fastq.gz`
do
fname=${fpath%_1.fastq.gz}
trimmomatic PE -threads 4 \
${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:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:20
done
mv *_cleaned_*.gz ${PROJECT_PATH}/cleaned_fastq/QC 過程が正しく行われていることを確認するために、クリーニング後の FASTQ に対して、もう一度 QC レポートを作成して確認するとよい。
マッピング
次にクリーニングしたリードをリファレンス配列にマッピングする必要がある。リードが RNA-Seq 由来ならば STAR、DNA-Seq 由来ならば BWA を使用することが推奨されている。ここで使用するサンプルは DNA-Seq 由来のリードであるため、BWA を使用する。
まず、BWA の使い方にしたがって、リファレンス配列のインデックスを作成する。
cd ${PROJECT_PATH}/genome
bwa index -p TAIR10_bwa Arabidopsis_thaliana.TAIR10.dna.toplevel.fa次に、クリーニングした FASTQ それぞれに対してマッピングを行う。なお、後ほど GATK でクオリティ補正を行うために、マッピングする際にリードグループ（RG）情報を与える必要がある。このリードグループの行に ID、PL、および SM の 3 種類のデータを代入する。ID には BAM ファイルの固有 ID を入れる。PL にはシーケンサーの種類を入れる。なお、シーケンサーの種類を表すキーワードはすで定義されている。Illumina 社のシーケンサーであれば ILLUMINA を入力する。また、SM はサンプル名を入れる。1 つのサンプルを複数回に分けてシーケンシングした場合は、ID は異なるが、SM が同じになることもある。
cd ${PROJECT_PATH}/cleaned_fastq
for fpath in `ls *_cleaned_1.fastq.gz`
do
fname=${fpath%_cleaned_1.fastq.gz}
bamRG="@RG\tID:"${fname}"\tPL:ILLUMINA\tSM:"${fname}
bwa mem -t 2 -R ${bamRG} ${PROJECT_PATH}/genome/TAIR10_bwa \
${fname}_cleaned_1.fastq.gz ${fname}_cleaned_2.fastq.gz > ${fname}.sam
done
mv *.sam ${PROJECT_PATH}/bam/BWA でマッピングした結果は、SAM フォーマットで保存される。これをコンピューターが効率よく扱えるように、リードを並べ替えてバイナリーファイルである BAM フォーマットに変換する。
cd ${PROJECT_PATH}/bam
for fpath in `ls *.sam`
do
samtools sort -@ 4 ${fpath} > ${fpath%.sam}.bam
doneバリアントコール
マッピング後、GATK を使用して、BAM ファイルからバリアントを検出する。この variant calling のプロセスは、前処理（duplication markup, BQSR）、SNPs/indels calling、そしてフィルタリングの 3 ステップからなる。これらの処理を一つずつ進めていく。なお、GATK 3.x のとき、前処理のプロセスとして local realignment も存在していたが、GATK4 ではこの処理がハプロタイプ推定時に行われるように実装されているため、前処理として実行する必要がない。
まず先に、リファレンス配列に対して、GATK 用のインデックスを作成する。また、リファレンス配列をこれから GATK の必須オプションとして頻繁に使っていくことになるので、ここで変数として宣言しておく。
cd ${PROJECT_PATH}/genome
REF=${PROJECT_PATH}/genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa
samtools faidx ${REF}
gatk CreateSequenceDictionary -R ${REF} -O ${REF%.fa}.dict
# gatk --java-options "-Xmx4G" CreateSequenceDictionary -R ${REF} -O ${REF%.fa}.dictなお、gatk コマンドを実行する時に、実行が途中でとまり、エラーが出力されないような場合は、メモリー不足である可能性がある。その際に、--java-options で Java が使えるメモリーを少し増やすと実行できるようになる。
前処理（duplication markup）
シーケンサーでシーケンシングしたリードの中に、PCR duplication などの重複リードを含む場合がある。このような重複リードは、variant calling の結果に含まれる偽陽性を高めてしまうことがある。そこで、variant calling に先立ち、これらの重複リードに目印をつけておく。目印をつけておけば、GATK はそのような重複リードを無視できるようになる。次のコードは、すべての BAM ファイルに対して重複リードのマークアップを行い、その結果を新しい BAM (*_markdup.bam) に保存するようにしてある。
cd ${PROJECT_PATH}/bam
for fpath in `ls *.bam`
do
gatk MarkDuplicates -I ${fpath} \
-O ${fpath%.bam}_markdup.bam \
-M ${fpath%.bam}_markdup_metrics.txt
samtools index ${fpath%.bam}_markdup.bam
done必要であれば、BAM ファイル中に含まれるアラインメントの統計量を計算することもできる。
cd ${PROJECT_PATH}/bam
for fpath in `ls *.bam | grep -v _markdup`
do
fname=${fpath%.bam}
gatk CollectAlignmentSummaryMetrics -I ${fname}_markdup.bam \
-R ${REF} \
-O ${fname}_alignment_metrics.txt
gatk CollectInsertSizeMetrics -I ${fname}_markdup.bam \
-O ${fname}_insert_metrics.txt \
-H ${fname}_insert_size_histogram.pdf
samtools depth -a ${fname}_markdup.bam > ${fname}_depth.txt
done前処理（BQSR）
BAM ファイル（あるいは FASTQ ファイル）中の記載された各リードのクオリティスコアは、シーケンサーから出力されたスコアである。このスコアを出力メカニズムはブラックボックスとなっているうえ、PCR エラーなどが反映されていない。そのため、variant calling のときに、間違ったクオリティで variant を検出してしまわないように、BAM ファイル中のクオリティを補正しておく必要がある。この補正作業を Base Quality Score Recalibratio (BQSR) とよぶ。
BQSR は、既知の variants を除いて残った部分のアラインメント領域のデータを利用して、機械学習モデルを構築し、そのモデルを使って全体のクオリティスコアを補正する。そのため、BQSR を行うには、既知の variants 情報を必要とする。既知の variants 情報を入手できる場合はその情報を使用すればよいが、入手できない場合は自ら暫定な variants 情報を作成する必要がある。
既知 variants 情報がある場合の BQSR
ヒトやマウスなどの生物では、既知の variants 情報がデータベースで保存されたり、あるいは先行研究としてまとめられている場合がある。それらの情報を vcf ファイルとして入手できれば、BQSR を次のように行うことができる。
gatk BaseRecalibrator -R ${REF} -I input.bam \
--known-sites dbsnp_146.b38.vcf.gz \
--known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
--known-sites resources_broad_hg38_v0_1000G.phase3.integrated.sites_only.hg38.vcf \
-O ${fname}_recal_data.table
gatk ApplyBQSR -R ${REF} -I input.bam -bqsr ${fname}_recal_data.table -O input_bqsr.bam既知 variants 情報がない場合の BQSR
既知の variants 情報がない場合は、クオリティ補正していない BAM データから一度 variant calling を行い、ここで得られる信用性の低い variants を既知情報として用いて BQSR を行う。では、サンプルデータを使用して、既知 variants 情報がない場合の BQSR を行う方法を示していく。
まず、それぞれの個体（BAM ファイル）に対してハプロタイプの推定を行う。ここでは gVCF ファイルが得られる。このファイルは VCF ファイルと同じフォーマットであるが、variants が存在しない領域の情報も含まれており、情報量が大きい。
cd ${PROJECT_PATH}/bam
for fpath in `ls *_markdup.bam`
do
fname=${fpath%_markdup.bam}
gatk HaplotypeCaller -R ${REF} --emit-ref-confidence GVCF \
-I ${fname}_markdup.bam -O ${fname}.g.vcf
done
mv *.g.vcf ../bqsr/
mv *.g.vcf.idx ../bqsr/次に、各個体の推定ハプロタイプを一つのファイルにまとめて、ジェノタイピングを行う（joint genotyping）。この処理で gVCF ファイルから variants 情報抽出され、VCF ファイルが作られる。この VCF ファイルには SNPs や indels などの情報が含まれる。
cd ${PROJECT_PATH}/bqsr
# list up all gVCF files
gvcf_files=""
for gvcf_file in `ls *.g.vcf`
do
gvcf_files=${gvcf_files}"-V ${gvcf_file} "
done
# build an database and merge gVCF files
echo -e "1\n2\n3\n4\n5" > intervals.list
gatk GenomicsDBImport -R ${REF} ${gvcf_files} -L intervals.list \
--genomicsdb-workspace-path gvcfs_db
gatk GenotypeGVCFs -R ${REF} -V gendb://gvcfs_db -O merged.vcf次に、ジェノタイピングにより得られた VCF ファイルの中から SNPs の情報だけを取り出して、別のファイルに保存し、フィルタリングを行う。これにより、信頼度のやや高い SNPs 情報を得ることができる。
cd ${PROJECT_PATH}/bqsr
gatk SelectVariants -R ${REF} -V merged.vcf --select-type-to-include SNP -O merged_snps.vcf
gatk VariantFiltration -R ${REF} -V merged_snps.vcf -O merged_snps_filtered.vcf \
-filter "QD < 2.0" --filter-name "QD2" \
-filter "QUAL < 30.0" --filter-name "QUAL30" \
-filter "SOR > 4.0" --filter-name "SOR4" \
-filter "FS > 60.0" --filter-name "FS60" \
-filter "MQ < 40.0" --filter-name "MQ40" \
-filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
-filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8"次に、indels に対しても同様な作業を行う。SNPs と indels を異なる閾値でフィルタリングする必要があるので、それぞれ分けて処理を行なった。
cd ${PROJECT_PATH}/bqsr
gatk SelectVariants -R ${REF} -V merged.vcf --select-type-to-include INDEL -O merged_indels.vcf
gatk VariantFiltration -R ${REF} -V merged_indels.vcf -O merged_indels_filtered.vcf \
-filter "QD < 2.0" --filter-name "QD2" \
-filter "QUAL < 30.0" --filter-name "QUAL30" \
-filter "FS > 200.0" --filter-name "FS200" \
-filter "SOR > 10.0" -filter-name "SOR10" \
-filter "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20"これにより、正確とは言えないが、暫定な SNPs 情報 (merged_indels_filtered.vcf) および indels 情報 (merged_indels_filtered) が得られる。この 2 つの情報を事前情報として、それぞれの BAM ファイルに対して BQSR を行う。
cd ${PROJECT_PATH}/bam
for fpath in `ls *_markdup.bam`
do
fname=${fpath%_markdup.bam}
gatk BaseRecalibrator -R ${REF} -I ${fname}_markdup.bam \
--known-sites ${PROJECT_PATH}/bqsr/merged_snps_filtered.vcf \
--known-sites ${PROJECT_PATH}/bqsr/merged_indels_filtered.vcf \
-O ${fname}_recal_data.table
gatk ApplyBQSR -R ${REF} -I ${fname}_markdup.bam -bqsr ${fname}_recal_data.table -O ${fname}_bqsr.bam
doneここまでの作業で、クオリティ補正後の *_bqsr.bam ファイルが得られる。なお、ここでは、BQSR を一回しか実行しなかったが、結果が収束するまでこの流れを繰り返して BQSR を実行することで、より信頼性の高い variant を得ることができる。
BQSR の効果確認
BQSR 実行前と実行後に、BAM ファイル中のクオリティがどのように改善されたかを調べる場合は、次のようにコードを実行する。ここで得られる結果は、あとで使用することはないので、BQSR の効果を確認という意味で実行すると良い。。
cd ${PROJECT_PATH}/bam
for fpath in `ls *_bqsr.bam`
do
fname=${fpath%_bqsr.bam}
gatk BaseRecalibrator -R ${REF} -I ${fname}_bqsr.bam \
--known-sites ${PROJECT_PATH}/bqsr/merged_snps_filtered.vcf \
--known-sites ${PROJECT_PATH}/bqsr/merged_indels_filtered.vcf \
-O ${fname}_recal_data.table.2
gatk AnalyzeCovariates -before ${fname}_recal_data.table -after ${fname}_recal_data.table.2 \
-plots ${fname}_recalibration_plots.pdf
doneバリアントコール
BQSR 後の BAM ファイル（*_bqsr.bam）に対して、variant calling を行なっていく。これ以降の操作は、BQSR 過程において暫定的な variant calling を行うときの操作と全く同じものとなっている。まず、各個体（BAM ファイル）に対してハプロタイプの推定を行う。
cd ${PROJECT_PATH}/bam
for fpath in `ls *_bqsr.bam`
do
fname=${fpath%_bqsr.bam}
gatk HaplotypeCaller -R ${REF} --emit-ref-confidence GVCF \
-I ${fname}_bqsr.bam -O ${fname}.g.vcf
done
mv *.g.vcf ../vcf/
mv *.g.vcf.idx ../vcf/次に、各個体の推定ハプロタイプをマージして、joint genotyping を行う。この処理によって得られる merge.vcf には SNPs や indels などが含まれている。また、それらの variants のクオリティは様々である。必要に応じてフィルタリングを行い、信頼性の高い variants を抽出して解析に持っていく。
cd ${PROJECT_PATH}/vcf
# list up all gVCF files
gvcf_files=""
for gvcf_file in `ls *.g.vcf`
do
gvcf_files=${gvcf_files}"-V ${gvcf_file} "
done
echo -e "1\n2\n3\n4\n5" > intervals.list
gatk GenomicsDBImport -R ${REF} ${gvcf_files} -L intervals.list \
--genomicsdb-workspace-path gvcfs_db
gatk GenotypeGVCFs -R ${REF} -V gendb://gvcfs_db -O merged.vcfフィルタリング
ジェノタイピング結果から SNPs の情報だけを抽出しフィルタリングを行う。
cd ${PROJECT_PATH}/vcf
gatk SelectVariants -R ${REF} -V merged.vcf --select-type-to-include SNP -O merged_snps.vcf
gatk VariantFiltration -R ${REF} -V merged_snps.vcf -O merged_snps_filtered.vcf \
-filter "QD < 2.0" --filter-name "QD2" \
-filter "QUAL < 30.0" --filter-name "QUAL30" \
-filter "SOR > 4.0" --filter-name "SOR4" \
-filter "FS > 60.0" --filter-name "FS60" \
-filter "MQ < 40.0" --filter-name "MQ40" \
-filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
-filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8"indels についても同様に抽出してフィルタリングを行う。
cd ${PROJECT_PATH}/vcf
gatk SelectVariants -R ${REF} -V merged.vcf --select-type-to-include INDEL -O merged_indels.vcf
gatk VariantFiltration -R ${REF} -V merged_indels.vcf -O merged_indels_filtered.vcf \
-filter "QD < 2.0" --filter-name "QD2" \
-filter "QUAL < 30.0" --filter-name "QUAL30" \
-filter "FS > 200.0" --filter-name "FS200" \
-filter "SOR > 10.0" -filter-name "SOR10" \
-filter "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20"なお、ここでは SNPs/indels のフィルタリングを行う際に、特定の統計量に対して閾値を設けてフィルタリングを行なった（hard-filtering）。このほかに、データベースあるいは先行研究から SNPs/indels に関する正しいと保証された事前知識を入手できる場合は、その事前知識で機械学習モデルを訓練して、そのモデルを使用してフィルタリングする方法もある。機械学習を利用するフィルタリングは Variant Quality Score Recalibration (VQSR) などと呼ばれている。