BWA 的基本介绍和主要算法
BWA 的基本介绍和主要算法
1. BWA 简介
BWA（Burrows-Wheeler Aligner）是一款将 DNA 序列比对到参考基因组上的工具，特别适用于低差异性的序列。它被广泛应用于基因组学和生物信息学研究，特别是在高通量测序数据处理和全外显子组测序（WES）分析中。
用途(比对)
BWA 可以将 DNA 序列映射到大型参考基因组
2. BWA 的三种算法
- BWA-backtrack：适用于**短于 100 bp **的 Illumina 测序读段。主要用于早期版本的 Illumina 数据。
- BWA-SW：适用于较长的序列（从 70 bp 到数百万 bp）。这种算法能够很好地处理较长的拼接片段或大型序列。
BWA-MEM：适用于序列长度从** 70 bp 到 1 Mbp** 的比对，最新且最常用的算法，适用于处理长读段和复杂序列。它通过寻找最大扩展匹配（MEM）来定位种子，然后使用 Smith-Waterman 算法进行局部比对，支持 indels（插入和缺失）处理。
BWA-MEM 特点：
- 能够处理长度从 70 bp 到 1 Mbp 的读段。
- 通过寻找 MEM 来进行种子化，支持 indels。
- 使用 Smith-Waterman 算法来扩展比对，精度和效率都很高。
3. BWA 的索引构建
在使用 BWA 进行比对之前，必须对参考基因组（如人类基因组）进行索引构建。BWA 提供了两种不同的索引构建算法，分别是 IS 和 bwtsw。参考基因组数据通常以 .fa 或 .fasta 格式提供。
基本命令
bwa index [-p prefix] [-a algoType] <in.db.fasta>
-p prefix：为生成的索引文件指定输出文件的前缀（可选）。-a algoType：指定构建索引的算法类型。可选的值有is或bwtsw。<in.db.fasta>：输入的参考基因组文件（FASTA 格式）。
算法选择
-
IS算法（IS linear-time algorithm）：- 是一种线性时间算法，用于构建后缀数组（suffix array）。
- 内存需求：需要参考基因组大小的 5.37 倍内存。
- 限制：不能处理大于 2 GB 的数据库。因此，
IS算法适用于小型参考基因组，但不适用于像人类基因组这样的大型基因组。 - 优点：简单且速度适中，是 BWA 的默认算法。
-
bwtsw算法（BWT-SW Algorithm）：- 适用于处理大型基因组数据，如整个人类基因组。
- 没有
2 GB的大小限制，因此更适合大规模基因组的数据处理。
生成的索引文件
索引构建完成后，BWA 会生成以下 5 个文件：
.bwt：Burrows-Wheeler Transform 结果。.pac：压缩后的参考基因组序列。.ann：注释信息。.amb：二进制的 N 区域标记文件（表示不确定碱基的位置）。.sa：后缀数组。
bwa index genome.fasta
genome.fasta.amb
genome.fasta.ann
genome.fasta.bwt
genome.fasta.pac
genome.fasta.sa
这些文件将用于后续的序列比对。
BWA 的比对
4. BWA 的序列比对
在完成参考基因组的索引构建后，可以使用 BWA 进行序列比对。根据不同的序列长度和数据类型，可以选择不同的算法来进行比对。在实际应用中，BWA-MEM 是最常用的比对算法，适用于大多数测序数据。
BWA-MEM 算法
BWA-MEM 适用于长度在** 70 bp 到 1 Mbp** 之间的序列，并且能处理长读段和复杂序列。BWA-MEM 通过寻找最大扩展匹配（MEM）来进行种子化，然后使用 Smith-Waterman 算法（https://blog.csdn.net/2302_80012625/article/details/142904050?spm=1001.2014.3001.5501）扩展这些种子，并支持 indels 的处理。
基本用法
bwa mem [options] <idxbase> <in1.fq> [in2.fq] > output.sam
[options]：一系列可选参数。<idxbase>：参考基因组的索引文件前缀。<in1.fq>和[in2.fq]：输入的质控后的 FASTQ 文件，in1.fq是必选的（单端测序或双端测序的第一个读段），in2.fq是可选的（双端测序的第二个读段）。output.sam：输出文件的路径，结果会以 SAM 格式输出。
示例
bwa mem -t 4 reference.fa sample_1.fq.gz sample_2.fq.gz > output.sam
# 使用 BWA 进行比对，并生成 SAM 文件
bwa mem -t 8 \ # 使用 8 个线程并行处理
-R '@RG\tID:H1\tSM:H1\tPL:illumina' \ # 添加 read group 信息：ID 为 H1，样本名为 H1，测序平台为 illumina
../01.ref/genome.fasta \ # 参考基因组文件（FASTA 格式） （已经通过 `bwa index` 构建）
/public/home/2022099/yiyaoran/workspace/02.ReadsFilter/H1_1.fq.gz \ # 第一条测序读取（FASTQ 格式，压缩）
/public/home/2022099/yiyaoran/workspace/02.ReadsFilter/H1_2.fq.gz \ # 第二条测序读取（FASTQ 格式，压缩）
-o H1.sam \ # 输出结果到 H1.sam 文件中
2>H1.bwa.log && # 将 BWA 比对过程中产生的日志输出到 H1.bwa.log 文件中
# 使用 samtools fixmate 处理配对比对信息
samtools fixmate -m H1.sam H1.fixmate.bam && # 将 H1.sam 修正后输出为 BAM 格式文件，并确保 mate 信息正确 (-m 选项)
# 对 BAM 文件进行排序
samtools sort -@ 2 -m 2G \ # 使用 2 个线程，分配每个线程 2G 内存进行排序
-o H1.fixmate.sort.bam \ # 将排序后的结果输出为 H1.fixmate.sort.bam
H1.fixmate.bam && # 输入修正后的 BAM 文件进行排序
# 使用 samtools markdup 标记重复比对
samtools markdup H1.fixmate.sort.bam H1.fixmate.sort.markdup.bam # 标记重复比对的片段，并将结果保存为 H1.fixmate.sort.markdup.bam
常用选项
-t：指定线程数，用于并行处理，加快比对速度。一般建议设定为 4。-k：设定最小种子长度。默认值是 19，较小的值会提高比对灵敏度，但降低速度。-w：指定比对的带宽，影响比对算法的搜索范围。默认值是 100。-R：为输出的 SAM 文件添加读取组（Read Group）信息，特别适合后续使用 GATK 进行分析。
Read Group 的重要性
在 BWA 的比对过程中，使用 -R 选项可以为输出的 SAM 文件添加 Read Group 信息。Read Group 信息包含了样本、测序平台、文库等重要信息，后续在使用 GATK 等工具进行分析时非常关键。
Read Group 示例
-R "@RG\tID:group1\tSM:sample1\tLB:lib1\tPL:ILLUMINA"
ID：Read Group 的 ID，一般为测序的 lane ID。SM：样本 ID。LB：文库名称。PL：测序平台（如ILLUMINA）。
调整灵敏度和速度的常用选项：
参考https://blog.csdn.net/2302_80012625/article/details/142904660?sharetype=blogdetail&sharerId=142904660&sharerefer=PC&sharesource=2302_80012625&spm=1011.2480.3001.8118
-
-k：设置最小种子长度。默认值为 19，增加该值可以减少误配率，但也会减少比对的灵敏度。bwa mem -k 19 ref.fa reads.fq > aln.sam -
-w：Band width，默认值为 100，设置为更大的值可以找到更长的插入或缺失片段。 -
-a：输出所有可能的比对，而不仅仅是最佳比对。bwa mem -a ref.fa reads.fq > aln.sam
比对得分调整：
参考:https://blog.csdn.net/2302_80012625/article/details/142904784
你可以调整 BWA-MEM 的比对得分参数来控制比对的准确性。
bwa mem -A 1 -B 4 -O 6 -E 1 ref.fa reads.fq > aln.sam
-A：匹配得分，默认值为 1。-B：错配罚分，序列错误率大约为{0.75 * exp[-log(4) * B/A]}，默认值为 4。-O：空位开启罚分，默认值为 6。-E：间隙延伸罚分，默认值为 1。