【直播】我的基因组 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

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-12

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏QQ音乐前端团队专栏

自己动手打造前端性能监控系统

我们从三个各方面,前端上报,数据收集和入库,数据展示来介绍了如何打造一个测速系统。

57710
来自专栏程序员的诗和远方

2018-0701_ARTS_week01

982
来自专栏企鹅号快讯

浅谈反馈式按钮的设计与实现

原创声明 ? 前言 前一段时间在网上闲逛看一些交互案例,偶然的看到几篇关于反馈式交互设计的文章,其中强调了反馈式设计的分类、重要性和机制,让我觉得在目前所负责的...

1827
来自专栏更流畅、简洁的软件开发方式

【角色】——分离开代码和权限需求,即实现代码和权限需求的解耦。

今天突然来了一个灵感,记录一下。以前总觉得说不清楚,看看这种表达方式是否可以说清。 两个原则:依赖接口编程,不要依赖实现编程;最小获知原则。 面向对象最重要的是...

2135
来自专栏hbbliyong

opoa介绍

一 定义       One Page, One Application(后面缩写为OPOA,或者1P1A), 含义很简单:一个页面就是一个应用。不再使用ifr...

2647
来自专栏向治洪

Android性能优化

讲到Android开发,就不得不谈一下Android的优化,不管是平时开发中我们需要注意的一些Android对Java的一些类的优化,还是实际开发中对性能的优化...

2076
来自专栏程序员的诗和远方

20180701_ARTS_week01

Two Sum Given an array of integers, return indices of the two numbers such that ...

641
来自专栏java学习

Java每日一练(2017/6/1)

Java基础 | 数据库 | Android | 学习视频 | 学习资料下载 课前导读 ●回复"每日一练"获取以前的题目! ●答案公布时间:为每期发布题目的第二...

2497
来自专栏Java帮帮-微信公众号-技术文章全总结

HTML5学习-day01【悟空教程】

网页超文本应用技术工作小组是一个以推动网络HTML 5 标准为目的而成立的组织。在2004年,由Opera、Mozilla基金会和苹果这些浏览器厂商组成。

773
来自专栏吉浦迅科技

PGI OpenACC 2018版:原来你是这样的编译器

对于CUDA Fortran用户来说,PGI编译器是必然要用到的。 其实PGI编译器不仅仅可以支持Fortran,还可以支持C/C++。而对于集群用户来说,要将...

4317

扫码关注云+社区