前面我们扯到bam文件的各种操作,vcf文件的各种操作,基础知识不牢固的同学可能已经云里雾里了。这次我们来讲一个简单的。就是拿到了fastq的测序数据,如何把全基因组分析给跑一遍。(不谈细节!) 首先就是fastq文件比对到参考基因组变成sam文件: head -40 read1.fq >tmp/read1.fq 一个简单的管道即可,如果管道不能确认是对的,就像我上面那样先拿一个小本文文件测试一下。由下图可以看到我们sort的bam文件不是按照染色体的1,2,3排序,而是按照chr10,chr11,,,,chr1,,chr2这样的顺序,这个对很多其它软件会不友好。 不过没关系,我们不跑GATK,这个bam文件足够了! 事实上,对我们真实的WGS数据来说,这一步耗时很严重的!(时间开销在后面) 第二个步骤,就是call variation咯,下面两个软件都可以,用起来也很简单。 : samtools mpileup -ugf ~/reference/genome/hg19/hg19.fa jmzeng.sorted.bam |bcftools call -vmO z -o jmzeng.bcftools.vcf.gzhead -40 read1.fq >tmp/read1.fq
|
|