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

相关文章

来自专栏Ceph对象存储方案

Luminous版本PG 分布调优

Luminous版本开始新增的balancer模块在PG分布优化方面效果非常明显,操作也非常简便,强烈推荐各位在集群上线之前进行这一操作,能够极大的提升整个集群...

3265
来自专栏张善友的专栏

Mix 10 上的asp.net mvc 2的相关Session

Beyond File | New Company: From Cheesy Sample to Social Platform Scott Hansel...

2627
来自专栏Golang语言社区

【Golang语言社区】GO1.9 map并发安全测试

var m sync.Map //全局 func maintest() { // 第一个 YongHuomap := make(map[st...

4858
来自专栏大内老A

The .NET of Tomorrow

Ed Charbeneau(http://developer.telerik.com/featured/the-net-of-tomorrow/) Exciti...

32710
来自专栏一个会写诗的程序员的博客

Spring Reactor 项目核心库Reactor Core

Non-Blocking Reactive Streams Foundation for the JVM both implementing a Reactiv...

2232
来自专栏跟着阿笨一起玩NET

c#实现打印功能

2932
来自专栏张善友的专栏

Miguel de Icaza 细说 Mix 07大会上的Silverlight和DLR

Mono之父Miguel de Icaza 详细报道微软Mix 07大会上的Silverlight和DLR ,上面还谈到了Mono and Silverligh...

2737
来自专栏陈仁松博客

ASP.NET Core 'Microsoft.Win32.Registry' 错误修复

今天在发布Asp.net Core应用到Azure的时候出现错误InvalidOperationException: Cannot find compilati...

4878
来自专栏张善友的专栏

LINQ via C# 系列文章

LINQ via C# Recently I am giving a series of talk on LINQ. the name “LINQ via C...

2675
来自专栏java 成神之路

使用 NIO 实现 echo 服务器

4827

扫码关注云+社区