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

相关文章

来自专栏CSDN技术头条

MongoDB开发版本3.1.8发布

MongoDB 3.1.8版本已发布。值得注意的是此次3.1.8作为开发版本,并不适用于生产环境中使用。接来下的3.2系列版本将供广大用户作为生产环境中使用,敬...

1996
来自专栏mini188

技术笔记:Indy的TIdSMTP改造,解决发送Html和主题截断问题

使用Indy来发邮件坑不少啊,只不过有比没有好吧,使用delphi6这种老工具没办法,只能使用了新一点的Indy版本9,公司限制。。。 1、邮件包含TIdTex...

1766
来自专栏FreeBuf

Meterpreter免杀技巧分享(亲测有效)

0×01 meterpreter简介 MetasploitFramework是一个缓冲区溢出测试使用的辅助工具,也可以说是一个漏洞利用和测试平台,它集成了各种平...

28610
来自专栏FreeBuf

基于unicorn-engine的虚拟机的实现(WxSpectre)

反病毒虚拟机是一个很有优势的工具,可以说反病毒软件是否存在模拟器是衡量反病毒软件能力的一个指标。反病毒虚拟机不光是内嵌在反病毒软件内部,来动态执行样本。这种虚拟...

3156
来自专栏性能与架构

使用Redis统计活跃用户

统计活跃用户这个案例非常经典,也是我当时学习redis时,接触到的第一个让我眼睛一亮的使用方式 场景 用户登录后需要记录,以便以后进行登录统计 统计需求主...

2846
来自专栏瓜大三哥

伪路径

1.什么是伪路径? 存在于设计之中,之所以叫伪路径,是因为这样的路径并未发挥真正的功能,在时序分析时不需要我们去分析。 2.为什么设置伪路径? 去除无效的时序...

2138
来自专栏点滴积累

Fish Shell

今天看到阮一峰同学的一篇博客(Fish shell 入门教程),讲述的非常详细、清楚,有兴趣的可以直接转去查看此文,本文仅提供一下个人使用心得。 一、fish ...

3196
来自专栏xiaoheike

es suggest did you mean资料

term suggester 根据提供的文档提供搜索关键词的建议,也就是关键词自动纠错。该链接介绍如何使用 term suggester 语法。term sug...

502
来自专栏杨建荣的学习笔记

一则数据库无法重启的案例分析(r8笔记第96天)

今天一个开发的同事找到我,说有个问题想咨询一下我,突然想起他昨天让我帮他处理一个工单,他这么一问我才想起来还没做,结果他说是另外一件事,说有个开 发测试的环境,...

3336
来自专栏菩提树下的杨过

Spring官网下载dist.zip的几种方法

Spring官网改版后,很多项目的完整zip包下载链接已经隐掉了,虽然Spring旨在引导大家用更“高大上”的maven方式来管理所依赖的jar包,但是完全没想...

1697

扫描关注云+社区