序列比对之BWA

2023-11-22 16:15:02 浏览数 (2)

文章

文献标题：Fast and accurate short read alignment with Burrows-Wheeler transform 发表日期：2009年7月 期刊：Bioinformatics 作者&单位：李恒 & Sanger Institute（桑格研究所） DOI：https://doi.org/10.1093/bioinformatics/btp324

简要用途

BWA是一个用于将DNA序列（特别是低差异性序列）映射到大型参考基因组（例如人类基因组）上的工具。它在基因组学和生物信息学研究中尤为重要。因为它能有效处理高通量测序数据，常常集成于WES分析流程，被广泛应用于基因组学研究，如在寻找与疾病相关的基因变异、理解种系发育关系等领域。

BWA包括三种主要算法：BWA-backtrack、BWA-SW和BWA-MEM。

  • BWA-backtrack 适用于短至100bp的Illumina序列读取，它是针对短序列设计的。
  • BWA-SW 和 BWA-MEM 适用于更长的序列，从70bp到数百万bp。BWA-MEM是最新的算法，为长读取和更复杂的序列提供了更高的精度和效率。

目前在以二代测试为主流的情况下，一般比对选择 BWA-MEM 算法就好。

性能特性：

  • BWA在速度和准确性方面表现出色，特别是在处理大型基因组数据时。它支持插入和缺失（indels）的比对，这是许多其他工具不提供的。
  • 与其他工具（如Bowtie）相比，BWA在处理更复杂的比对情况（如允许indels）时可能稍慢，但在精确性方面更为可靠。

如何安装

一般来说我们默认推荐是使用conda来管理软件

代码语言：javascript复制
conda create -n wes ##先创建小环境，如果已经创建，无需再次执行
conda activate wes
codna install -y bwa

二进制安装

当然，采用二进制安转也很简单

代码语言：javascript复制
wget -c https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.17.tar.bz2
tar -xvf bwa-0.7.17.tar.bz2 -C /home/yjzhang/biosoft
cd ~/biosoft/bwa-0.7.17/
make
./bwa mem 
ln -s ~/biosoft/bwa-0.7.17/bwa ~/biosoft/mybin/

最小化使用

bwa 软件的作用是将序列比对到参考基因组上，在比对之前，首先需要对参考基因组建立索引，再比对

bwa 的使用需要两种输入文件：

  • Reference genome data（fasta格式 .fa, .fasta, .fna）
  • Short reads data (fastq格式 .fastaq, .fq)

1构建索引

代码语言：javascript复制
bwa index [-p prefix] [-a algoType] <in.db.fasta>

# -p 输出文件前缀 [可选参数]
# -a 构建索引采用的算法「可选参数」;可选择指定采用 is 还是 bwt-sw 算法构建索引

##example
bwa index -a bwt-sw reference.fa

关于构建索引算法选择

  1. is (IS linear-time algorithm)：
    • 这是一种用于构建后缀数组的线性时间算法。
    • 它需要的内存量是数据库大小的5.37倍，其中N代表数据库的大小。
    • IS算法的速度适中，但有一个重要的限制：它不能处理大于2GB的数据库。这意味着对于非常大的基因组数据，如整个人类基因组，这种方法可能不适用。
    • 尽管有这个限制，IS算法由于其简单性，被设定为默认算法。
  2. bwtsw (Algorithm implemented in BWT-SW)：
    • 与IS算法不同，bwtsw算法适用于那些需要处理大规模基因组数据的场景，比如处理整个人类基因组这样的大型基因组。

索引建立好之后，会生成5个文件，后缀分别为

  • bwt
  • pac
  • ann
  • amb
  • sa

2比对

代码语言：javascript复制
## 基本用法
~/biosoft/bwa-0.7.17/bwa mem  -t 4 /public/reference/index/bwa/hg38 
/public/project/father/sampel_1.fq.gz 
/public/project/father/sample_2.fq.gz 
-o ~/out-bwa/sample.sam

# -t 线程，默认为1，增加线程数，会加速程序的运行时间，但也会增加内存和CPU资源的使用，因此也并不是设置的越大越好，以我们经验来说建议设置为4个线程即可，如果设置超过8个线程，对于运行速度的提升来说基本是微乎其微；
# -o 设定输出路径及文件
# 默认情况下，输入一个序列文件认为是单端测序，输入两个序列文件则是双端测序

其余用法及参数

采用不同算法，比对命令也会不同

BWA-MEM 算法

BWA-MEM算法适用于序列长度在70bp到1Mbp之间的序列比对。算法首先通过寻找MEMs来种子化（seeding）比对。MEMs是指在参考基因组中能找到的与查询序列完全匹配的最长片段。这些MEMs作为潜在比对位置的初始点。找到MEMs后，BWA-MEM算法使用带有affine-gap惩罚的Smith-Waterman算法来扩展这些种子。Affine-gap惩罚是一种在序列比对中用于处理插入和缺失（indels）的技术。Smith-Waterman算法是一种经典的动态规划算法，用于局部序列比对，能够找到最优的局部比对。

其对应的子命令为mem, 用法如下：

代码语言：javascript复制
Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]
[options] #是一系列可选的参数。
<idxbase> #参考基因组的索引文件，上面通过bwa index构建好的文件便是。
<in1.fq>和[in2.fq] #输入的是质控后的fastq文件。

## 调整灵敏度和速度
bwa mem -k 19 -w 100 ref.fa reads.fq > aln.sam
-k ##设定最小种子长度。长度小于此值的匹配将被忽略。这个值影响比对的灵敏度和速度。默认值19
-w ##Band width。大于此值的间隙将无法被找到。不过最大间隙长度还受到评分矩阵和击中长度的影响，并不是完全由此项决定。默认值100

bwa mem -a ref.fa reads.fq > aln.sam
-a ## 参数使得所有可能的比对结果都会输出，而不仅是最佳比对。

##调整比对打分参数
bwa mem -A 1 -B 4 -O 6 -E 1 ref.fa reads.fq > aln.sam
-A ##匹配得分，默认值1。
-B ##错配罚分。序列错误率大约为：{.75 * exp[-log(4) * B/A]}。
-O ##空位罚分。
-E ##间隙延伸罚分。长度为k的间隙的得分为O   k*E（即-O是为了开启一个零长度的间隙）。

## 实例应用
sample='7E5239'
bwa mem -t 5 -R "@RGtID:$sampletSM:$sampletLB:WGStPL:Illumina" /public/biosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38 /home/wes/clean/7E5239.L1_1_val_1.fq.gz /home/yjzhang/wes/clean/7E5239.L1_2_val_2.fq.gz -o ~/wes/align

-t ##设定线程数，默认为1
-R #Read Group的字符串信息，这是一个非常重要的信息，以@RG开头，它是用来将比对的read进行分组的。不同的组之间测序过程被认为是相互独立的，这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。在Read Group中，有如下几个信息非常重要：
- 1) ID，这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；
- 2) PL，指的是所用的测序平台，这个信息不要随便写！特别是当我们需要使用GATK进行后续分析的时候，更是如此！这是一个很多新手都容易忽视的一个地方，在GATK中，PL只允许被设置为：ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这几个信息。基本上就是目前市场上存在着的测序平台，当然，如果实在不知道，那么必须设置为UNKNOWN，名字方面不区分大小写。如果你在分析的时候这里没设置正确，那么在后续使用GATK过程中可能会出错
#我们上面的例子用的是PL:illumina。如果你的数据是CG测序的那么记得不要写成CG！而要写COMPLETE。
- 3) SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本；
- 4) LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB；
除了以上这四个之外，还可以自定义添加其他的信息，不过如无特殊的需要，对于序列比对而言，这4个就足够了。这些信息设置好之后，在RG字符串中要用制表符（t）将它们分开。

BWA-backtrack 算法

对应的子命令为aln/samse/sample

bwa aln

寻找输入读取的后缀数组（SA）坐标，并在此过程中设定了允许的差异（即不匹配）的数量。

  • SA坐标（Suffix Array coordinates）：
    • 后缀数组是一种数据结构，用于存储字符串的所有后缀的排序列表。在序列比对中，它被用于快速定位（或“映射”）输入序列（即读取）在较长参考序列（如人类基因组）中的位置。
    • 寻找输入读取的SA坐标意味着确定这些读取在参考基因组中可能对应的位置。
  • 最大差异数（maxSeedDiff和maxDiff）：
    • 这部分说明了在确定SA坐标时允许的最大差异数，差异通常指的是读取和参考序列之间的不匹配。
    • maxSeedDiff：这是在读取的第一个子序列（或“种子”）中允许的最大差异数。这个子序列的长度由seedLen参数确定。这意味着在进行初步的比对（种子比对）时，序列间允许有一定数量的不匹配。
    • maxDiff：这是在整个读取序列中允许的最大差异数。这意味着在整个读取和参考序列的比对中，允许的不匹配总数不应超过这个数值。
代码语言：javascript复制
## 单端比对
bwa aln ref.fa reads.fq > reads.sai

## 双端比对
bwa aln ref.fa reads1.fq > reads1.sai
bwa aln ref.fa reads2.fq > reads2.sai

## 其余用法
bwa aln -t 4 ref.fa reads.fq > reads.sai
-t ## 设置线程

bwa aln -l 32 ref.fa reads.fq > reads.sai
-l ## 调整比对时使用的种子长度。这对于处理特定类型的测序数据或优化比对性能可能很有用

bwa aln -q 20 ref.fa reads.fq > reads.sai
-q ## 设定一个质量阈值，该命令会从每个读取的3'端开始，去除质量分数低于20的碱基，然后对修剪后的读取进行比对。这个参数可以提高比对准确性，但是过度修剪可能导致有效数据的损失，因此需要根据实际数据质量和分析目标谨慎设置这个参数。

bwa aln -n 0.03 -o 1 ref.fa reads.fq > reads.sai
# 调整错配和间隙罚分
-n NUM：## 设置最大编辑距离（即允许的最大差异数）。如果是整数，则直接指定最大编辑距离；如果是浮点数，则表示在2%均匀碱基错误率下缺失比对的比例。对于不同长度的读取，最大编辑距离会自动选择。
-o INT ## 最大间隙开启次数。

## 使用BAM格式作为输入文件：
## 单端
bwa aln ref.fa -b reads.bam > reads.sai
## 双端
bwa aln ref.fa -b1 reads.bam > 1.sai
bwa aln ref.fa -b2 reads.bam > 2.sai

bwa samse

将单端测序数据的比对结果（存储在.sai文件中）和原始的测序读取（通常是.fastq格式）转换为标准的SAM（Sequence Alignment/Map）格式

代码语言：javascript复制
## 单端比对
bwa aln ref.fa reads.fq > reads.sai
bwa samse ref.fa reads.sai reads.fq > reads.sam

## 其余用法
bwa samse -n 5 ref.fa reads.sai reads.fq > reads.sam
-n ## 参数指定每个读取的最大额外比对数。这里设置为5，意味着每个读取最多记录5个额外比对位置。

bwa samse -r '@RGtID:group1tSM:sample1' ref.fa reads.sai reads.fq > reads.sam
-r ## 参数为SAM文件添加读取组信息

bwa sampe

将双端测序数据的比对结果和原始的测序读取转换为标准的SAM格式，在这个过程中，如果读取对在参考基因组中有多个可能的比对位置（即重复读取对），这些读取对将被随机分配到可能的比对位置之一

代码语言：javascript复制
## 基本用法
bwa aln ref.fa reads1.fq > reads1.sai
bwa aln ref.fa reads2.fq > reads2.sai
bwa sampe ref.fa reads1.sai reads2.sai reads1.fq reads2.fq > aln.sam

## 其余用法
bwa sampe -a 500 ref.fa reads1.sai reads2.sai reads1.fq reads2.fq > aln.sam
-a # 参数可以指定读取对之间的最大插入大小。这对于筛选出异常长度的成对读取很有用。

bwa sampe -o 2 ref.fa reads1.sai reads2.sai reads1.fq reads2.fq > aln.sam
-o ## 参数可以控制重复读取对的处理方式。这里的值2意味着重复读取对将被随机放置。

bwa sampe -n 10 -N 20 ref.fa reads1.sai reads2.sai reads1.fq reads2.fq > aln.sam
-n ## 限制正确配对的读取对在XA标签中输出的最大比对数，
-N ## 限制不正确配对的读取对的最大比对数

bwa sampe -r '@RGtID:group1tSM:sample1' ref.fa reads1.sai reads2.sai reads1.fq reads2.fq > aln.sam
-r ## 添加read 组信息

BWA-SW 算法

在in.fq文件中比对查询序列。当mate.fq文件存在时，执行双端（paired-end）比对。双端比对模式仅适用于Illumina短插入片段文库的读取。在双端比对模式下，BWA-SW可能仍会输出split alignments，但它们都会被标记为未正确配对；如果配对读取在基因组中有多个局部比对位置，其配对位置将不会被记录。

对应的子命令为bwasw, 基本用法如下 ：

代码语言：javascript复制
##单端比对
bwa bwasw ref.fa reads.fq > aln.sam

## 双端比对
bwa bwasw ref.fa reads1.fq reads2.fq > aln.sam

## 其余用法
bwa bwasw -t 4 ref.fa reads.fq > aln.sam
-t ## 设置线程

bwa bwasw -z 2 -w 100 ref.fa reads.fq > aln.sam
-z ## 选项可以调整比对的准确性（代价是速度），默认值是1 
-w ## 选项可以调整Band width，影响比对算法的搜索范围，默认值是33。

bwa bwasw -q 20 ref.fa reads.fq > aln.sam
-q ##指定最大间隙长度，这可以影响间隙的识别

bwa bwasw -T 50 ref.fa reads.fq > aln.sam
-T ## 设定最小得分阈值，只有达到这个阈值的比对才会被记录

多个算法该如何选择

不同算法适用不同数据，官方建议是：

  • 对于 Illumina、454和IonTorrent平台生成的单端读取（长度通常超过70bp）以及长度可达几兆碱基的组装拼接片段
代码语言：javascript复制
bwa mem ref.fa reads.fq > aln.sam
  • 对于长度短于 ~70bp 的Illumina单端测序读取：
代码语言：javascript复制
 bwa aln ref.fa reads.fq > reads.sai
 bwa samse ref.fa reads.sai reads.fq > aln-se.sam
  • 对于 Illumina、454和IonTorrent平台生成的双端读取（长度通常超过70bp）
代码语言：javascript复制
bwa mem ref.fa read1.fq read2.fq > aln-pe.sam
  • 对于长度短于 ~70bp 的Illumina双端测序读取：
代码语言：javascript复制
bwa aln ref.fa read1.fq > read1.sai
bwa aln ref.fa read2.fq > read2.sai
bwa sampe ref.fa read1.sai read2.sai read1.fq read2.fq > aln-pe.sam
  • 对于 PacBio 或Oxford Nanopore 测试读取比对到参考基因组：
代码语言：javascript复制
bwa mem -x pacbio ref.fa reads.fq > aln.sam
bwa mem -x ont2d ref.fa reads.fq > aln.sam

为什么bwasw 目前很少使用？

bwasw 专门设计用于处理长读取数据。它适用于以下情境：

  1. 长读取数据：bwasw效果最佳的是对于长度超过100bp的长读取，尤其是来自454、Ion Torrent和早期Illumina平台的读取。
  2. 组装拼接片段：对于数千到数百万碱基对长度的组装拼接片段，bwasw也是一个合适的选择。

但是随着测序技术的发展，尤其是随着第三代测序平台（如PacBio和Oxford Nanopore）的兴起，测序读取的长度大大增加。这些新平台产生的长读取对比对算法提出了新的挑战，例如更高的错误率和更复杂的错误模式。因此，针对这些长读取优化的比对算法（如Minimap2）变得更为流行，它们能更有效地处理长读取数据的特点。而对于第二代测序技术，如Illumina平台，它们通常产生较短的读取（75-300bp），并且具有较低的错误率。对于这类读取，bwa mem算法表现更优，因为它在处理短读取方面更为高效和准确。所以目前bwa bwasw在常规流程中的使用很少。

参考

  • bwa使用手册：https://github.com/lh3/bwa
ref 数据 算法 索引 线程

0 人点赞