如何用 samtools 和 bcftools call snp

samtools 之前博文已經介紹過一些常用的方法。本篇主要說下如何利用samtools 和 bcftools來call snp。

和其他工具一樣，bam檔案都要經過處理（另見博文）。假如對C17樣本進行call snp, 資料為：

LC17-1_L002.sorted.rmp.rg.recal.bam

LC17-2_L006.sorted.rmp.rg.recal.bam

LC17-3_L002.sorted.rmp.rg.recal.bam

RC17-1_L003.sorted.rmp.rg.recal.bam

RC17-2_L004.sorted.rmp.rg.recal.bam

RC17-3_L004.sorted.rmp.rg.recal.bam

資料準備好後，執行指令：

samtools mpileup -P ILLUMINA -f RAP_cDAN.fasta -EgD \

LC17-1_L002.sorted.rmp.rg.recal.bam LC17-2_L006.sorted.rmp.rg.recal.bam \

LC17-3_L002.sorted.rmp.rg.recal.bam RC17-1_L003.sorted.rmp.rg.recal.bam \

RC17-2_L004.sorted.rmp.rg.recal.bam RC17-3_L004.sorted.rmp.rg.recal.bam \

> samtools_result.bcf

指令解釋：

mpileup 是samtools中call snp的工具。

-P 指platform, 現在短reads測序一般是ILLUMINA。

-f 後跟參考序列，序列檔案必須提前建好index。

-E， Extended BAQ（base alignment quality） computation,

如果有的話，會提高檢測出MNPs的靈敏度，當然會輕微的減下特異度。

-g Compute genotype likelihoods and output them in the binary call

format(BCF).

-D Output per-sample read depth

> 是将結果儲存到samtools_result.bcf檔案中

最終得到的samtools_result.bcf 是二進制檔案，到此完成了call snp的第一步。

得到bcf檔案以後，第二步執行指令：

bcftools view -cNegv samtools_result.bcf > samtools_result.vcf

veiw 是bcftools中主要的方法，‘Convert between BCF and VCF ， call variant candidates

and estimate allele frequencies.’

-c  Call variants using Bayesian inference .

-N Skip sites where the REF field is not A/T/G/C.

一般的參考基因組序列都是由四種堿基組成，不知道還能有什麼？ 難道是沒出來部分的N 麼 ？ 也可以不加這個參數，

   我測試過，這種情況非常非常少。

-e 其實也可以不加， 因為如果前面有-c 那麼-e就預設被使用了。‘Perform max-likelihood inference

only,including estimating the site allele frequency, testing Hardy-Weinberg

equlibrium and testing associations with LRT’.

-g Call per-sample genotypes at variant sites.(用-c的方法)

-v output variant site only .我們當然隻關心編譯位點

> 将結果儲存到samtools_result.vcf檔案中。 是vcf格式的，具體關于vcf檔案格式解釋見其他博文。