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

相关文章

来自专栏程序猿DD

我最常用的Intellij IDEA快捷键

原文:http://www.radcortez.com/my-most-useful-intellij-idea-keyboard-shortcuts/ 译文:...

1918
来自专栏james大数据架构

Eclipse快捷键大全

一、实用类快捷键 1 常用熟悉的快捷键 CTRL+C(复制)、CTRL+X(剪切)、CTRL+Z(撤销)、CTRL+F(查找)、CTRL+H(搜索文件或字符串)...

17310
来自专栏守望轩

Visual Studio 2008 每日提示(三十六)

# 361、按Ctrl+K, Ctrl+v在解决方案快速查找对象 原文链接: Ctrl+K, Ctrl+v allows you to quickly sea...

3497
来自专栏康怀帅的专栏

CoreOS 容器 Rkt 简单介绍

由于 Docker 已经成为事实上的容器老大,这里暂且将 rkt 内容放入 docker 文件夹。哈哈 官方网站:https://coreos.com/rkt/...

3987
来自专栏腾讯云实验室

腾讯云实验室的正确投稿姿势

腾讯云实验室上线了在线投稿能力,现在除了官方推出的实验室之外,允许有能力的开发者把自己的技术和经验通过在实验室投稿的方式来进行分享和传播。

3962
来自专栏Ldpe2G的个人博客

Mxnet Scala Package 学习笔记 一

882
来自专栏开源FPGA

FPGA计算3行同列数据之和

实验:FPGA计算3行同列数据之和 实验要求:PC机通过串口发送3行数据(一行有56个数据,3行共有56*3=168个数据)给FPGA,FPGA计算3行同一列数...

1888
来自专栏沈唁志

Composer进阶使用之版本约束表达式的使用

例如我们想要下载5.1版本的ThinkPHP包,我们可以通过composer.json文件:

372
来自专栏Lambda

Idea 常用快捷键

———–自动代码——– 常用的有fori/sout/psvm+Tab即可生成循环、System.out、main方法等boilerplate样板代码 例...

1916
来自专栏Ldpe2G的个人博客

Mxnet Scala Package 学习笔记 一

从刚开始接触Mxnet这个框架到现在已经大概两年了。MXNet最吸引我的地方就是它提供了

813

扫码关注云+社区