samtools(bcftools)を使ってvariant call

ヒト以外のモデル生物とかモデル生物じゃないもののvariant callingをしたいときは、samtools mpileupを使ってvcfを作りましょう。

準備

samtools
bcftools
vcfutils.pl

アプリのダウンロードはhttp://www.htslib.org/から
インストールは上記ページに書いてある通りに
(make -> ./configure -> make install)
vcfutils.plは追記

reference(下の例ではref.fa)はsamtools faidxでインデックスをつけておく
(いらないかも）

samtools faidx ref.fa

コマンド

Samtools mpileup -> bcftools call でvariantを抽出

SNPと短いINDELが取れてくるはず

variant_call.sh

samtools mpileup -uf <ref.fa> -t AD,INFO/AD,DP,DV,DPR,INFO/DPR,DP,DP4,SP -v <hoge1.bam> <hoge2.bam> | bcftools call -c -v -o hoge_variant.vcf.txt

vcfutils.pl varFilter hoge_variant.vcf.txt > filterd_hoge_variant.vcf.txt

samtools mpileupのオプション

-u ;uncompressed 出力
-f file ;リファレンスファイル指定
-t list ;追加するタグを指定
-v ;VCFを出力（このオプション無しでもbcftoolsで読めるようだけど試したことはない）

bcftools callのオプション

-c ;original calling method (multi allelic, rare variantは-mを使うように、と説明には書いてある）
-v ;VCFを出力
-o ;出力先ファイル名指定

samtoolsのマニュアル
"*.vcf"だけだと住所録ファイルのvcf形式として認識されてしまうので、後ろに　".txt"をつけておく

bcftoolsでmpileupもできるらしい

bcftools mpileup -Ou -f <ref.fa> <sample1.bam> <sample2.bam> <sample3.bam> | bcftools call -vmO z -o <study.vcf.gz>

bcftools mpileupのオプション

-Ou ;uncompressed BCF出力
-f file ;リファレンスファイル指定

positionごとのGenotypeカウント

filterd_hoge_variant.vcf.txtをRで読みこむ

multi sample callingの時、判定されたgenotypeのカウントをします

0/0 -> refと同じ
0/1 -> ヘテロ
1/1 -> ホモ

準備(必要なパッケージ)

VariantAnnotation
Bioconductorより

Rのコード

genotype_count.r

library(VariantAnnotation)

in_f <- "filterd_hoge.vcf.txt"
out_f <- "count_gt_filterd_hoge.txt"

vcf <- readVcf(in_f,"hoge")　

GT <- geno(vcf)$GT　
ct <- apply(GT,1,table)
none <- sapply(ct,"[", "0/0")
hetero <- sapply(ct,"[","0/1")
homo <- sapply(ct,"[","1/1")
x <- cbind(row.names(GT),none,hetero,homo)
write.table(x,file=out_f,quote=F,sep="\t",row.names=F,col.names=T)