|
在这里，我们可以考虑比较对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： 我-O选择的是z格式输出，没有添加多线程，所以对43G的bam文件处理了30个小时。 -v 就是说，我只输出有变异的位点，而不是全基因组坐标的输出。 同时也可以选择多线程计算 因为全基因组的计算量很大,可以指定 --threads 24 便是使用24线程 。在先跑通最简单命令的基础之上，不要一成不变，可以通过读软件的文档来理解参数的意义，通过适当的调整参数，来更好地解决你的计算任务，但是参数不要随便改，因为不够robust的算法在面对不同复杂度的数据的时候，很有可能会带给你错误的结果 参考： https://samtools./bcftools/bcftools.html http://cbsu.tc./lab/doc/Variant_workshop_Part1.pdf 文：Jimmy、阿尔的太阳 图文编辑：吃瓜群众
|
|
|
