在Window7平台利用samtools尋找SNPs / INDELs @ 有勁的基因資訊

在先前的部落格中我們介紹過了在Windows 7平台下利用bowtie或bowtie2將NGS定序的資料和reference序列作alignment，並利用samtools進行檔案的處理以及利用IGV觀看瀏覽alignment結果

詳細內容可參照：

今天我們再來教大家如何在Windows 7下利用samtools來找出定序資料中SNP (single-nucleotide polymorphism) 以及INDEL (insertion or deletion)。

前置作業：準備好以下檔案

1. 安裝samtools (可參照bowtie2教學中步驟五部分)
2. 定序資料跟reference的alignment檔案 (bam format, sorted)
3. Reference序列fasta檔案，並以samtools faidx建立index
(指令: samtools faidx reference.fa，此步驟會產生reference.fa.fai檔案)

如下圖中，本次demo使用的reference序列為lamda_virus.fa，alignment檔案為lambda.sorted.bam

20131204_pic1

步驟一：利用samtools mpilup指令以建立每個base的堆疊資料
在Windows的命令提示字元視窗中輸入以下指令:
samtools mpileup –gf [reference.fasta] [aln.sorted.bam] > mpileup.bcf

如下圖所示，本教學所使用的reference為lambda_virus.fa；alignment檔案為lambda.sorted.bam，此兩個檔案均存放於mpileup_demo資料夾下。另外output的結果我們將它寫在mpileup.bcf檔案裡

20131204_pic2

步驟二：，並利用bcftools view找出SNP/INDEL
由於前一步驟產生的檔案僅僅是根據alignment檔案建立每個base的堆疊資料，我們還需要利用samtools裡的另外一支程式”bcftools”來找出SNP/INDEL。bcftools 指令如下：
bcftools view –vcg mpileup.bcf > var.raw.vcf
(如下圖)

20131204_pic3

註1：此步驟產生的檔案可直接在命令提示字元視窗下用more指令觀看或是以NotePad++、Excel等程式開啟。

註2：如對於vcf檔案格式各欄位代表意義有疑問的話可參考以下表格

20131204_pic4

20131204_pic5

步驟三：利用vcfutils.pl進一步過濾SN/INDEL結果
如果你的電腦已經有安裝Perl的話，可以利用samtools裡的一支程式vcfutils.pl來過濾前一步驟所得到的SNP/INDEL結果。例如：如果我們要對於SNP calling所需要的sequencing depth設上限為100條reads的話，可在命令提示字元視窗中輸入如下指令：
perl vcfutils.pl varFilter –D100 var.raw.vcf > var.flt.vcf

20131204_pic6