【直播】我的基因组 43：简单粗糙的WGS数据分析流程

前面我们扯到bam文件的各种操作，vcf文件的各种操作，基础知识不牢固的同学可能已经云里雾里了。这次我们来讲一个简单的。就是拿到了fastq的测序数据，如何把全基因组分析给跑一遍。(不谈细节！)

首先就是fastq文件比对到参考基因组变成sam文件：

head -40 read1.fq >tmp/read1.fq 
head -40 read2.fq >tmp/read2.fq
~/biosoft/bwa/bwa-0.7.15/bwa  mem -t 20 -M ~/reference/index/bwa/hg19 read1.fq read2.fq| samtools sort -O bam -o jmzeng.sorted.bam

一个简单的管道即可，如果管道不能确认是对的，就像我上面那样先拿一个小本文文件测试一下。由下图可以看到我们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 
~/biosoft/freebayes/freebayes/bin/freebayes  -f  ~/reference/genome/hg19/hg19.fa  jmzeng.sorted.bam >jmzeng.freebayes.vcfhead -40 read2.fq >tmp/read2.fq

大家非常后台留言最多的就是这3个步骤的耗时问题

分别是84290、134143，76406 秒！也就是23.4hours、37.3hours、21.2hours，正好一个双休日，哈哈，完美！

两个call variation的步骤是并行的。也就是说完成一个全基因组数据(300G的原始数据)的分析，是需要整整两天两夜的！

但是大家可能在朋友圈多次看到各种宣传贴21小时完成千人的全基因组分析，为什么呢？是因为硬件条件的不同，他们有着相当大的计算资源。他们的内存和存储空间都要比我们自己所用的计算资源大不知道多少倍。同时，还有各样本并行处理及GPU的加速，所以运行速度那么快也就不足为奇了。

文：Jimmy

图文编辑：吃瓜群众