搜索
查看: 2623|回复: 1

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-7 19:39:35 | 显示全部楼层 |阅读模式
【直播】我的基因组 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

图文编辑:吃瓜群众






本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:【直播】我的基因组 42: 不同lane的variation的比较
下一篇:【直播】我的基因组 44:比对文件画profile和heatmap图
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

0

主题

15

帖子

88

积分

注册会员

Rank: 2

积分
88
发表于 2017-2-19 01:13:58 | 显示全部楼层
之前贴已经建议过哈,可以用sambamba多线程节省samtools的大量时间。
另外freebayes的准确性和GATK calling的准确性如何(GATK里面best practice中各种realign, recalicate很耗时)?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 09:35 , Processed in 0.029199 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.