【直播】我的基因组70:比对文件并不能完美的还原出测序文件

前面我们说到过可以用软件或者自己写脚本从已经比对到参考基因组的sam/bam格式文件提取出原始的测序fastq文件。

但是我在IGV里面检查bam文件的时候发现了一些难以理解的现象,所以趁这个机会把它们探究清楚。

bwa工具的不同版本影响大吗? bwa对同样测序文件同样参数比对多次结果一样吗? bwa对一个PE reads只输出两条记录吗? bwa的-M参数是干什么的? bwa会截断原始fastq序列吗?

首先可以找到各个版本的软件的下载地址,如下:

https://sourceforge.net/projects/bio-bwa/files/?source=navbar

## 安装软件很简单,我就不多说了 wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.8.tar.bz2 tar jxvf bwa-0.7.8.tar.bz2 cd bwa-0.7.8 make

一般公司都会用已经建立好的流程,不太会频繁更新软件,所以是VN:0.7.8-r455 版本的bwa工具,而我喜欢下载最新版,是VN:0.7.15-r1140

我制作了一个PE的fastq测序数据,如下:

## 首先是read1.fq ST-E00142:330:H33KVALXX:5:1219:28280:57301 /1 ACGCCTGTAATCCCAGCACTTTGGGAGGCCAAGGCAGGCGGATCACCTGAGGTTAGGAGTTCAAGACCACCCTGGCCAACATGGTGAAACCCCATCTCTACTAAAAAACAAAAAATTAACTAGGTGTGGTGGCAGGCGCCTGTAATCCCA + AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJA-FJJAJJJJF<A-FJJA-<F<F-<FA7F77-F-77FA7F-7FJ<A---7AA7-F--AJ-FJJAAA<<AF-JJJFJJJJJJJFJJJJJJJ<JJJJFJJJFJJJJJJJJJJFJFFFJ ## 然后是read2.fq ST-E00142:330:H33KVALXX:5:1219:28280:57301 /2 CCAGGAGGCAGAGGGTGCCAGTGATCCGAGATTGCGCCATTGCACTCCTTCTTAATGAAACGGCGCCCACCGCGATCTACACACGAACCCTACGCGCCGCTCTTCCGATCTACGCCTGTAATCCCAGCACTTTGGGAGGCCAAGGTGGGC + A-JF<)7FJAF<A7))<)<-<<JAFFAJFF--JJJFF-F-JF77F77-----A7--77-777F-7--7J7A7-7FF-JJ7<------7777A7----7F-77-7--<---7J<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA

你们需要把两个fastq文件用上面的序列自己制作好;

然后测试两个版本的bwa结果,代码如下:

/home/jianmingzeng/biosoft/bwa/bwa-0.7.15/bwa mem -M /home/jianmingzeng/reference/index/bwa/hg19 read1.fq read2.fq 1>bwa-0.7.15.sam /home/jianmingzeng/biosoft/bwa/bwa-0.7.8/bwa mem -M /home/jianmingzeng/reference/index/bwa/hg19 read1.fq read2.fq 1>bwa-0.7.8.sam

下面的内容我虽然复制粘贴在这里,但是我建议你弄到notepad++等编辑器里面仔细观看,最重要的是,你自己走一般这个过程,不然你根本不知道我在说什么。

VN:0.7.15-r1140比对结果是:

@PG ID:bwa PN:bwa VN:0.7.15-r1140 CL:/home/jianmingzeng/biosoft/bwa/bwa-0.7.15/bwa mem -M /home/jianmingzeng/reference/index/bwa/hg19 read1.fq read2.fq

ST-E00142:330:H33KVALXX:5:1219:28280:57301 65 chr1 1346332 60 150M = 1346519 188 ACGCCTGTAATCCCAGCACTTTGGGAGGCCAAGGCAGGCGGATCACCTGAGGTTAGGAGTTCAAGACCACCCTGGCCAACATGGTGAAACCCCATCTCTACTAAAAAACAAAAAATTAACTAGGTGTGGTGGCAGGCGCCTGTAATCCCA AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJA-FJJAJJJJF<A-FJJA-<F<F-<FA7F77-F-77FA7F-7FJ<A---7AA7-F--AJ-FJJAAA<<AF-JJJFJJJJJJJFJJJJJJJ<JJJJFJJJFJJJJJJJJJJFJFFFJ NM:i:0 MD:Z:150 AS:i:150 XS:i:92

ST-E00142:330:H33KVALXX:5:1219:28280:57301 129 chr1 1346519 44 48M102S = 1346332 -188 CCAGGAGGCAGAGGGTGCCAGTGATCCGAGATTGCGCCATTGCACTCCTTCTTAATGAAACGGCGCCCACCGCGATCTACACACGAACCCTACGCGCCGCTCTTCCGATCTACGCCTGTAATCCCAGCACTTTGGGAGGCCAAGGTGGGC A-JF<)7FJAF<A7))<)<-<<JAFFAJFF--JJJFF-F-JF77F77-----A7--77-777F-7--7J7A7-7FF-JJ7<------7777A7----7F-77-7--<---7J<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA NM:i:1 MD:Z:14T33 AS:i:43 XS:i:35 SA:Z:chr1,120410765,-,42M108S,0,0; XA:Z:chr3,+57932288,17M1I30M102S,2;

ST-E00142:330:H33KVALXX:5:1219:28280:57301 401 chr1 120410765 0 42M108H = 1346332 -119064475 GCCCACCTTGGCCTCCCAAAGTGCTGGGATTACAGGCGTAGA AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ<J7-- NM:i:0 MD:Z:42 AS:i:42 XS:i:42 SA:Z:chr1,1346519,+,48M102S,44,1;

VN:0.7.8-r455 比对结果是:

@PG ID:bwa PN:bwa VN:0.7.8-r455 CL:/home/jianmingzeng/biosoft/bwa/bwa-0.7.8/bwa mem -M /home/jianmingzeng/reference/index/bwa/hg19 read1.fq read2.fq

ST-E00142:330:H33KVALXX:5:1219:28280:57301 65 chr1 1346332 60 150M = 1346519 188 ACGCCTGTAATCCCAGCACTTTGGGAGGCCAAGGCAGGCGGATCACCTGAGGTTAGGAGTTCAAGACCACCCTGGCCAACATGGTGAAACCCCATCTCTACTAAAAAACAAAAAATTAACTAGGTGTGGTGGCAGGCGCCTGTAATCCCA AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJA-FJJAJJJJF<A-FJJA-<F<F-<FA7F77-F-77FA7F-7FJ<A---7AA7-F--AJ-FJJAAA<<AF-JJJFJJJJJJJFJJJJJJJ<JJJJFJJJFJJJJJJJJJJFJFFFJ NM:i:0 MD:Z:150 AS:i:150 XS:i:0

ST-E00142:330:H33KVALXX:5:1219:28280:57301 129 chr1 1346519 44 48M102S = 1346332 -188 CCAGGAGGCAGAGGGTGCCAGTGATCCGAGATTGCGCCATTGCACTCCTTCTTAATGAAACGGCGCCCACCGCGATCTACACACGAACCCTACGCGCCGCTCTTCCGATCTACGCCTGTAATCCCAGCACTTTGGGAGGCCAAGGTGGGC A-JF<)7FJAF<A7))<)<-<<JAFFAJFF--JJJFF-F-JF77F77-----A7--77-777F-7--7J7A7-7FF-JJ7<------7777A7----7F-77-7--<---7J<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA NM:i:1 MD:Z:14T33 AS:i:43 XS:i:35 SA:Z:chr1,120410765,-,42M108S,0,0;

ST-E00142:330:H33KVALXX:5:1219:28280:57301 401 chr1 120410765 0 42M108H = 1346332 -119064475 GCCCACCTTGGCCTCCCAAAGTGCTGGGATTACAGGCGTAGA AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ<J7-- NM:i:0 MD:Z:42 AS:i:42 XS:i:42 SA:Z:chr1,1346519,+,48M102S,44,1;

我们都知道sam文件的第2列是flag,最正常应该是99和147,但是65和129是还算正常的PE reads的比对情况,分别代表左右两端

但是 401 这个东西我以前没有留意,查了一下,代表的是 not primary alignment

这里又涉及到了那个-M参数的意义:

The BWA-MEM algorithm performs local alignment. It may produce multiple primary alignments for different part of a query sequence. This is a crucial feature for long sequences. However, some tools such as Picard’s markDuplicates does not work with split alignments. One may consider to use option -M to flag shorter split hits as secondary.

而且里面还涉及到了H/S这两种比对结果:

我们也知道比对结果里面的H和S分别代表hard和soft的clipping,但是知道里面的具体含义的可能不多,如果是H模式,那么fastq会被截断后比对,这样bam文件里面看到的fastq就是一个部分序列,所以就不可能从bam文件里面还原出fastq序列啦。如果是S的话,虽然被截断的序列也是比对不说,但是在bam里面仍然会出现完整的fastq序列。

我这里已经回到了最开始我提出来的5个问题,我知道一般人看不懂!

我也没指望你们看懂,只有跟着做,动手的人才能看懂,就这样了,O(∩_∩)O谢谢

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

原文发表时间:2017-04-25

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏沈唁志

如何简单计算PHP网站是否已经最高负载

1975
来自专栏pangguoming

AngularJS 中文资料+工具+库+Demo 大搜集

中文学习资料: 中文资料且成系统的就这么多,优酷上有个中文视频。 http://www.cnblogs.com/lcllao/archive/2012/10/1...

3676
来自专栏IMWeb前端团队

后台:nodejs 前台:vue 全栈开发 外卖平台系统

关于 一直考虑写一个功能齐全的完整Nodejs项目,但苦于没有找到合适的类型,而且后台系统无法直观的感受到,需要有一个前台项目配合,因此迟迟没有动笔。恰好前一段...

2980
来自专栏农夫安全

绕过某商业waf实现getshell

0x00 前言 目标站:http://www.xxxx.com/ Aspcms 拿后台就不说了,太多姿势了。 主要说下后台getshell(过云盾拦截) A...

3415
来自专栏程序猿DD

Spring Cloud构建微服务架构:分布式服务跟踪(抽样收集)【Dalston版】

通过 TraceID和 SpanID已经实现了对分布式系统中的请求跟踪,而这些记录的跟踪信息最终会被分析系统收集起来,并用来实现对分布式系统的监控和分析功能,比...

3406
来自专栏java一日一条

干货!如何正确使用Git Flow

我们已经从SVN 切换到Git很多年了,现在几乎所有的项目都在使用Github管理, 本篇文章讲一下为什么使用Git, 以及如何在团队中正确使用。

1174
来自专栏用户画像

白盒测试与黑盒测试的定义与区别

白盒测试方法按照程序内部的结构测试程序,检验程序中的每条通路是否都能按预定要求正确工作,而不顾它的功能。

572
来自专栏项勇

[Android Studio]之设置字体大小与样式

2865
来自专栏CSDN技术头条

详解 NoSQL 数据库的分布式算法

系统的可扩展性是推动NoSQL运动发展的的主要理由,包含了分布式系统协调,故障转移,资源管理和许多其他特性。这么讲使得NoSQL听起来像是一个大筐,什么都能塞进...

2209
来自专栏iOSDevLog

Google Colab免费GPU教程

现在,你可以开发深度学习与应用谷歌Colaboratory -on的免费特斯拉K80 GPU -使用Keras,Tensorflow和PyTorch。

3755

扫码关注云+社区