BWA使用详解

参考来源http://starsyi.github.io/2016/05/24/BWA-%E5%91%BD%E4%BB%A4%E8%AF%A6%E8%A7%A3/

一、BWA 简介

BWA,即Burrows-Wheeler-Alignment Tool。BWA 的学习主要来自参考网站使用手册;BWA 是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包。它由三个不同的算法:

BWA-backtrack: 是用来比对 Illumina 的序列的,reads 长度最长能到 100bp。-
BWA-SW: 用于比对 long-read ,支持的长度为 70bp-1Mbp;同时支持剪接性比对。
BWA-MEM: 推荐使用的算法,支持较长的read长度,同时支持剪接性比对(split alignments),但是BWA-MEM是更新的算法,也更快,更准确,且 BWA-MEM 对于 70bp-100bp 的 Illumina 数据来说,效果也更好些。

二、 BWA 使用

在运用这三种算法之前,需要先利用 BWA 的 index 命令,构建出参考基因组的 FM-index,而对与上述的三种不同的算法而言,又有不同的命令:

aln/samse/sampe ----> BWA-backtrack (samse 中的 se 是 single-end 的简写,而 sampe 中的 pe 是 paired-end 的简写)。
bwasw ----> BWA-SW
mem ----> BWA-MEM
3.1 建立索引 index

在进行 reads 的比对前,需要对 fasta 文件构建 FM-index。用法和参数如下:

index Usage:
      bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>
OPTIONS: 
      -p STR   输出数据库的前缀;【默认和输入的文件名一致,输出的数据库在其输入文件所在的文件夹,并以该文件名为前缀。】
      -a [is|bwtsw]   构建index的算法,有两个算法: is 是默认的算法,虽然相对较快,但是需要较大的内存,当构建的数据库大于
               2GB的时候就不能正常工作了。 bwtsw 对于短的参考序列式不工作的,必须要大于等于10MB, 但能用于较大的基因组数据,比如人的全基因组。

#根据reference genome data(e.g. ref.fa) 建立 Index File例子:
$ bwa index ref.fa -p genome###可以不加-p genome,这样建立索引都是以ref.fa为前缀
3.2 mem比对

该算法先使用 MEM(maximal exact matches) 进行 seeding alignments,再使用 SW(affine-gap Smith-Waterman) 算法进行 seeds 的延伸。BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配,导致 reads 有多个最佳匹配位点。这对 long reads 的比对时比较重要的结果。但是却会和 Picard 的 markDuplicates 程序不兼容。使用方法如下:

mem Usage: bwa mem [options] ref.fa reads.fq [mates.fq]

简单的来讲,mem 使用的 MEMs(maximal exact matches) 进行seedling alignments, 再使用 SW(affine-gap Smith-Waterman 算法)进行seeding extending.
mem 进行局部比对,因此,对于一条序列的不同区域可能会产生多种最优匹配结果, 这对于long reads 来说尤为重要。 有些软件如 Picard’s markDuplicates 跟 mem 的这种剪接性比对不兼容,在这种情况下,可以使用 –M 选项来将 shorter split hits 标记为次优。
特别说明
如果 mates.fq 缺省,且参数 –p 未设定,那么 reads.fq 被认为是 single-end;如果 mates.fq 存在,且参数 –p 未设定,那么 mem 命令会认为 read.fq 和 mates.fq 中的 i-th reads 组成一个read对 (a read pair),这个模式是常用的 paired-end mode。如果参数 –p 被设定,那么, mem 命令会认为 read.fq 中的 第 2i-th 和 第 (2i + 1)-th 的 reads 组成一个 read 对 (a read pair),这种方式也被成为交错式的(interleaved paired-end)。 在这种情况下,即使有 mates.fq,也会被忽略。
常用的参数如下

-t   INT 线程数,默认是1。
-M   将 shorter split hits 标记为次优,以兼容 Picard’s markDuplicates 软件。
-p   若无此参数:输入文件只有1个,则进行单端比对;若输入文件有2个,则作为paired reads进行比对。若加入此参数:则仅以第1个文件作为输入(输入的文件若有2个,则忽略之),该文件必须是read1.fq和read2.fa进行reads交叉的数据。
-R   STR 完整的read group的头部,可以用 '\t' 作为分隔符, 在输出的SAM文件中被解释为制表符TAB. read group 的ID,会被添加到输出文件的每一个read的头部。
-T   INT   当比对的分值比 INT 小时,不输出该比对结果,这个参数只影响输出的结果,不影响比对的过程。-a 将所有的比对结果都输出,包括 single-end 和 unpaired paired-end的 reads,但是这些比对的结果会被标记为次优。

#例子:
$ bwa mem ref.fa reads.fq > mem-se.sam
$ bwa mem ref.fa read1.fq read2.fq > mem-pe.sam
$ bwa mem -t 4 -M -R "\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\" ref.fa read1.fastq read2.fastq > mem-pe.sam 2> ./mem-pe.log
3.3 align/samse/sampe比对

用法如下:

#对于single-read
bwa aln [options] ref.fa read.fq > aln_sa.sai
bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam
#对于pair-reads:
bwa aln [options] ref.fa read1.fq > aln_sa1.sai
bwa aln [options] ref.fa read2.fq > aln_sa2.sai
bwa sampe [options] ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam
3.4 bwasw

对输入的第1个文件的所有序列进行比对。如果输如有 2 个文件,则进行 paired-end 比对,此模式仅对 Illumina 的 short-insert 数据进行比对。在 Paired-end 模式下,BWA-SW依然输出剪接性比对结果,但是这些结果会标记为 not properly paired; 同时如果有多个匹配位点,则不会写入 mate 的匹配位置。常用参数如下:

bwasw Usage: bwa bwasw [ options ] ref.fasta reads.fq [mate.fq] > aln.sam
-t INT 使用的线程数

例子:
$ bwa bwasw genome long_read.fq > aln.sam$ bwa bwasw genome read1.fq read2.fq > aln-pe.sam
3.5 backtrack

经典的 bwa 先使用 aln 命令将单独的 reads 比对到参考序列,再使用 samse 或 sampe 生成 sam 文件。常用例子:

$ bwa aln genome read1.fq > aln_sa1.sai$ bwa aln genome read2.fq > aln_sa2.sai
$ bwa samse genome aln_sa1.sai read1.fq > aln_se.sam
$ bwa sampe genome aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln_pe.sam
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 158,425评论 4 361
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,058评论 1 291
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 108,186评论 0 243
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,848评论 0 204
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,249评论 3 286
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,554评论 1 216
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,830评论 2 312
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,536评论 0 197
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,239评论 1 241
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,505评论 2 244
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,004评论 1 258
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,346评论 2 253
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,999评论 3 235
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,060评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,821评论 0 194
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,574评论 2 271
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,480评论 2 267

推荐阅读更多精彩内容