原标题:还在为微生物重测序变异检测发愁?samtools帮助你!
Samtools作为一款操作序列比对结果文件(SAM/BAM)的工具,能够灵活转换sam/bam,并且能基于参考序列和Sam/bam文件进行变异位点检测,并通过bcftools进行变异结果统计。本期将基于上期bwa软件的比对结果,利用samtools进行变异检测。
一、软件安装
下载地址http://www.htslib.org/download/,利用wget下载得到samtools-1.5.tar.bz2,bcftools-1.5.tar.bz2;
解压缩:tar -jxvf bcftools-1.5.tar.bz2
tar -jxvf samtools-1.5.tar.bz2
安装(详细的安装信息可查看bcftools-1.5(samtools-1.5)文件夹下的INSTALL):
cd bcftools-1.5 ( samtools-1.5)
./configure –prefix=/where/to/install#设置安装路径
Make
Make install
安装完成后,我们后面所用到的可执行程序都在上面设置的安装路径(/where/to/install)的bin文件夹下;
配置环境变量
a.若要临时修改环境变量,可直接在终端输入下面一行命令:
Export PATH=/where/to/install/bin:$PATH
b.要永久修改环境变量可将下面第一行添加到~/.bash_profile(针对当前用户)或者/etc/profile(针对所有用户)文件的末尾,再执行第二行命令即可:
Export PATH=$PATH: /where/to/install/bin;
source ~.bash_profile 或者source /etc/profile
现在samtools 和bcftools已经安装好了,下面我们就进行calling SNPs和Indels了。
二、使用流程
1. 输入文件
比对结果文件sample.sort..bam和参考基因组序列ref.fa(上一期已对其建立索引,此处不再建立索引);
2. Variant Calling主要流程
Samtools index sample.sort.bam#对bam文件建立索引,默认生成文件为bam文件加.bai,此处生成sample.sort.bam.bai;
Samtools rmdup sample.sort.bam sample.dedup.bam#去除PCR等实验过程中产生的多余duplications;
Samtools mpileup –q 20 –d 8000 –ugo sample.bcf –f ref.fa sample.dedup.bam#variant calls
参数说明:
-q, --min-MQ mapQ最小值,低于这个值则跳过;
-d, --max-depth 最大覆盖深度;
-g, --BCF 生成bcf格式文件;
-o, --output FILE 输出文件;
-f , --fasta-ref FILE参考序列(fasta format)的索引文件名;
其他常用参数:
-A --count-orphans 不丢弃异常配对reads;
-l --positions FILE 包含区域位点的位置列表文件或者bed区域文件;
-r --region REG 在指定区域产生pipeup,和-l一起使用;
-I --skip-indel 不检测INDEL;
-m --min-ireads 候选INDEL的最小间隔的reads数;
-v --VCF 生成VCF格式文件;
-u --uncompressed 产生未压缩的vcf/bcf文件;
Bcftools call –vmO z –o sample.vcf.gz sample.bcf#将bcf转化为vcf文件,bcf为vcf的二进制格式
参数说明:
-v --variants-only 仅输出变异位点,不加此参数,则输入所有位点信息;
-c --consensus-caller the original calling method;
-m -- multiallelic-caller alternative model for multiallelic and rare-variant calling,calling的方法,与-c不能同时使用;
-O --output-type 指定输出文件类型, 'b'表示压缩后BCF文件; 'u'表示未压缩的BCF文件; 'z'表示压缩的VCF文件; 'v' 表示未压缩的VCF文件;
三、vcf文件的过滤统计
bcftools filter -O z -o sample.filtered.vcf.gz -i 'QUAL>20 && DP>5' sample.vcf.gz#过滤质量值低于20并且深度小于5的变异位点
参数说明:
-O --output-type 指定输出文件类型, 'b'表示压缩后BCF文件; 'u'表示未压缩的BCF文件; 'z'表示压缩的VCF文件; 'v' 表示未压缩的VCF文件;
-o, --output FILE 输出文件;
-i --Include 筛选出满足exp条件的变异位点;
-e --exclude 剔除满足exp条件的变异位点,与-i相反;
--threads 线程数;
Note: 上面的过滤条件中可以加入DP4参数,过滤出高质量的变异位点
Vcf 文件格式:由#开头的注释部分(图1)和变异主体部分(图2,图3),注释部分主要包括软件版本,命令行,参考序列信息以及主体部分参数的注释,主体部分由以tab隔开的8列组成,依次为参考序列名,variant的位置,variant的ID,参考序列的碱基,Phred格式(Phred_scaled)的质量值,过滤结果和第8列variant 的详细信息(图3),采用tag=value形式,tag间用“;”隔开;
图1 VCF文件注释部分
图2 VCF文件变异主体部分前7列
图3 VCF文件变异主体部分详细信息和基因型列
bcftools stats sample.filtered.vcf.gz > sample.stat.file# 统计SNPs和Indels变异数目等,包括转换,颠换数目;
sample.stat.file部分内容如下图所示,number of SNPs为SNPs总数目,number of indels为发生indels的总数目,ts表示转换数目,tv表示颠换数目:
bcftools view –v snps –g hom sample.filtered.vcf.gz# 仅输出纯合SNPs,通过改变参数可以输出杂合SNPs,纯合/杂合Indels;
参数说明:
-h 仅输出头部注释文件;
-H 不输出头部注释文件;
-O --output-type 指定输出文件类型, 'b'表示压缩后BCF文件; 'u'表示未压缩的BCF文件; 'z'表示压缩的VCF文件; 'v' 表示未压缩的VCF文件;
-T, --threads 设置线程数目
-i --Include 筛选出满足exp条件的变异位点;
-e --exclude 剔除满足exp条件的变异位点,与-i相反;
-v --types [snps|indels|mps|other],仅输出指定类型的变异位点信息;
-V --exclude-type [snps|indels|mps|other],仅输出不包含指定类型的变异位点信息;
-g --genotype [hom|het|miss],仅输出指定基因型的位点。
四、参考
Samtools官方说明文档:http://www.htslib.org/doc/#manual-pages
Bcftools官方说明文档:http://samtools.github.io/bcftools/bcftools.html
参考文献:
Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup (2009)The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics, 25, 2078-9. [PMID:19505943]
Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011 Nov 1;27(21):2987-93. Epub 2011 Sep 8. [PMID:21903627]doi:10.1093/bioinformatics/btr509;
Danecek P., Schiffels S., Durbin R. Multiallelic calling model in bcftools (-m) [link]
Li H. Improving SNP discovery by base alignment quality. Bioinformatics. 2011 Apr 15;27(8):1157-8. doi: 10.1093/bioinformatics/btr076. Epub 2011 Feb 13. [PMID:21320865]
Durbin R. Segregation based metric for variant call QC [link]
Li H, Mathematical Notes on SAMtools Algorithms [link]
The Genomes Project Consortium. An integrated map of genetic variation from 1,092 human genomes[J]. Nature, 2012, 491(7422):56-65.doi:10.1038/nature11632
Wu Y, Jing R, Dong Y, et al. Functional annotation of sixty-five type-2 diabetes risk SNPs and its application in risk prediction[J]. Scientific Reports, 2017, 7.
供稿:协云基因微生物事业部 韩娜
想了解更多?
那就赶紧来关注我们
有问题请联系协和基因云
微信ID:
geneworks返回搜狐,查看更多
责任编辑: