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

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

MENU

SNP情報の取得[Samtools mpileup]

Samtools mpileupの使い方

samtools, bcftools, vcftools

BAMファイルから変異情報を記述したVCF (Variant Call Format)という形式のファイルをbuildします。全ゲノムスケールの様々な解析にはVCFファイルが前提となることが多いため、GWAS(ゲノムワイド関連解析)には不可欠なステップになります。ブログ主の使ったことのあるSNP callingツールとしてSamtools_mpileupとGATKという米国Broad研究所が開発したソフトウェアがあります。GATKの使い方は、公式HPのBest Practiceをよく読んで実践する必要があるため、またの機会にいくつかの記事に分けて細かく紹介します。今回は、比較的簡易に実行可能なsamtools mpileupの使い方を紹介します。

本ブログでは、GATKの使い方を詳細かつ具体的に紹介したGATK解説シリーズを連載しています
GATKはsamtools pileupよりも統合的なツールで1つです。習得は大変ですが、使えるようになるメリットは非常に大きいので、是非挑戦してみてください。
ブログ主はGATKに挑戦される方を、ものすごく応援しています。

eupatho-bioinfomatics.hatenablog.com


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

なお、本ブログの過去記事↓

eupatho-bioinfomatics.hatenablog.com

で作成したマラリア原虫のBAMファイルを例として用います。

1. 必要なツールをインストールする。

conda install -c bioconda samtools
conda install -c bioconda vcftools
conda install -c bioconda bcftools

2. samtools mpileupを実行する。
workingdirには作業するディレクトリーのPATHを入力し、REFにはリファレンスゲノムのPATHを指定する。
BAMには、前回用意したBAMファイルを格納したフォルダを指定している。
また、ID_list.txtは過去記事で作成したAccession numberのリストである。

workingdir=/path/to/workingdir
REF=/$workindir/3D7_genomes/PlasmoDB-52_Pfalciparum3D7_Genome.fasta
BAM=/$workingdir/sort.bam
for ID in `cat ID_list.txt`
do
cd $workingdir
samtools mpileup --max-depth 100 --min-BQ 15 -guI --fasta-ref $REF $BAM/${ID}.sort.bam |  \
bcftools call  -cv -o $workingdir/${ID}.var.raw.bcf
bcftools view $workingdir/${ID}.var.raw.bcf | vcfutils.pl varFilter -d 10 > $workingdir/${ID}.vcf
done

samtools mpileupのオプション
--max-depth: カバレッジが非常に高い領域を処理するために必要なメモリと時間を減らすことができます。このオプションにゼロを渡すと、この制限が事実上なくなります。
--min-BQ: 考慮するbae qualityの最低値を規定する。
--fasta-ref: リファレンス配列のFastaファイル。
bcftools callのオプション
-c: original calling method
-v: VCFを出力する。
vcfutils.pl varFilterのオプション
-d: 最少カバレッジのフィルタリングオプション。

3. VCFファイルをmergeする。

for ID in `cat ID_list.txt`
do
cd $workingdir
bcftools view $ID.vcf -Oz -o $ID.vcf.gz
bcftools index $ID.vcf.gz
done

vcf-merge ERR1081238.vcf.gz ERR1081239.vcf.gz ERR1081241.vcf.gz ERR1081242.vcf.gz ERR1081254.vcf.gz ERR1099214.vcf.gz ERR1081255.vcf.gz ERR1081257.vcf.gz ERR1099215.vcf.gz ERR1081261.vcf.gz ERR1081262.vcf.gz ERR1081263.vcf.gz ERR1081264.vcf.gz ERR1081265.vcf.gz ERR1081283.vcf.gz ERR1081284.vcf.gz ERR1081285.vcf.gz ERR1081287.vcf.gz ERR1106549.vcf.gz ERR1081237.vcf.gz | bgzip -c > merge.vcf.gz

出力されたmerge.vcf.gzを使って、個別の解析を行う。

今回はこれで終わりです。よければ他の記事のも見ていってください。 なお、今回の記事の作成にあたり、samtools 使い方 mpileup ( calling SNPs ) & annotationVCFやBCF を扱う bcftoolsを参考にさせていただきました。

eupatho-bioinfomatics.hatenablog.com