More than 1 year has passed since last update.


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 -> ホモ


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


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)