【直播】我的基因组 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 条评论
登录 后参与评论

相关文章

来自专栏代码GG之家

google 进入分屏后在横屏模式按home键界面错乱( 四)

google 进入分屏后在横屏模式按home键界面错乱( 四) 你确定你了解分屏的整个流程? ? 代码阅读,请到此处http://androidxref.com...

24180
来自专栏生信技能树

给学徒的ATAC-seq数据实战

查看文章发现数据上传到了GEO,是:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66581

90240
来自专栏码匠的流水账

聊聊partition的方式

一般来说,数据库的繁忙体现在:不同用户需要访问数据集中的不同部分,这种情况下,我们把数据的各个部分存放在不同的服务器/节点中,每个服务器/节点负责自身数据的读取...

17110
来自专栏安恒信息

LOCKY勒索者新花样:通过PDF投递

摘 要 最近安恒APT团队截获一个新版的LOCKY勒索者病毒样本,区别之前大多数样本采用WORD文档投递并用宏代码远程下载执行的方式,该样本在原有的WORD文档...

30060
来自专栏生信技能树

基因组重测序的unmapped reads assembly探究 【直播】我的基因组86

在前面的直播基因组系列,我们讲解过那些比对不少我们人类的参考基因组序列的数据,其实可以细致的进行探究。 直播】我的基因组(十五):提取未比对的测序数据 这里主...

433160
来自专栏北京马哥教育

用python爬虫抓站的一些技巧总结

这些脚本有一个共性,都是和web相关的,总要用到获取链接的一些方法,再加上simplecd这 个半爬虫半网站的项目,累积不少爬虫抓站的经验,在此总结一下,那么以...

29150
来自专栏生信技能树

从WGS测序得到的VCF文件里面提取位于外显子区域的【直播】我的基因组84

首先要下载并且得到人类基因组的外显子坐标记录文件 这里我用的参考基因组版本仍然是hg19,所以去CCDS数据库里面下载对应版本,并且格式化成BED文件。 wge...

48690
来自专栏开发 & 算法杂谈

构建动态数据竞争检测平台

之前的文章分别介绍了基于Lockset算法、基于happens-before关系以及基于两者混合的方法。对于这些方法,已有的一些论文中提到的有关实验对比可能都比...

24440
来自专栏性能与架构

NoSql数据库的主要模型

KVP键值对模型 是一组两个关联的数据项,非常简单,有很高的灵活性和可扩展性 随着数据量的增加,KVP的计算也自然增加,所以使用KVP模型的数据库是指数型的 典...

33140
来自专栏机器学习算法与Python学习

值得收臧 | 从零开始搭建带GPU加速的深度学习环境(操作系统、驱动和各种机器学习库)

关键字全网搜索最新排名 【机器学习算法】:排名第一 【机器学习】:排名第一 【Python】:排名第三 【算法】:排名第四 原文:https://medium....

40160

扫码关注云+社区

领取腾讯云代金券