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で記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
なお、本ブログの過去記事↓
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 ) & annotation、VCFやBCF を扱う bcftoolsを参考にさせていただきました。