查看原文
其他

【直播】我的基因组25:用bcftools来call variation

2016-12-21 生信菜鸟团 生信技能树

在这里,我们可以考虑比较对3种不同的bam文件来分别call variation的差异,探索对bam文件的不同过滤模式对snp calling的影响,分别是,原始比对的bam文件,去除低质量和多比对还有PCR重复的bam,以及用GATK进行realign的bam文件。

首先采取最经典的软件来做这个工作,后面我们会使用其它软件并且比较:

代码如下:


samtools mpileup -ugf ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta \

P_jmzeng.filter.rmdup.bam  |bcftools call -vmO z -o P_jmzeng.rmdup.bcftools.vcf.gz


在这里  \ 是换行符 在脚本里面很多时候一条命令很长  包含数据的绝对路径以及命令的选项  比如GATK 的很多命令就很长 加入换行符还是一行命令 在脚本中看起来更加清晰


所以只需要理解samtools mpileup bcftools call的参数的意义就可以了


首先我们来看samtools的mpileup这个命令我选择的参数

很明显,我就选择了3个参数,-f 来输入有索引文件的fasta参考序列; -g 输出到bcf格式给bcftools软件来做input,而-u就是说不要压缩,因为我们要通过管道直接把数据给bcftools的。这个mpileup以前为pileup,如果大家看到的教程比较陈旧,会有这样的疑问。


接下来我们看bcftools:

同时也可以选择多线程计算  因为全基因组的计算量很大,可以指定 --threads 24 便是使用24线程 。在先跑通最简单命令的基础之上,不要一成不变,可以通过读软件的文档来理解参数的意义,通过适当的调整参数,来更好地解决你的计算任务,但是参数不要随便改,因为不够robust的算法在面对不同复杂度的数据的时候,很有可能会带给你错误的结果

参考:







文:Jimmy、阿尔的太阳

图文编辑:吃瓜群众




您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存