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

詳細內容可參照:

1. bowtie教學

2. bowtie2, samtools教學

3. IGV教學

今天我們再來教大家如何在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  


如要查看vcfutils.pl varFilter功能的其他參數預設值以及設定方法的話,可以直接在命令提示字元視窗中打入:perl vcfutils.pl varFilter直接按Enter。(如下圖)

 20131204_pic7  

 

整個samtools calling SNP/INDEL的步驟如上所示,只要簡單的三個步驟就可以達到目的!

註:目前windows 版本的samtools僅更新到0.1.12a,其關於指令的串接以及一些指令的output仍有bug存在,如果你按照此步驟操作有問題的話,請改用Linux或Mac電腦並安裝最新版本之samtools (目前最新版本為0.1.19)

 

 

 

Yourgene Bioscience  

 

YourGene 發表在 痞客邦 PIXNET 留言(0) 人氣()