526互联

生信：一起读官方文档 BWA 篇

发布时间 2023-08-29 17:48:48作者: 仗剑天涯横刀笑

一起读官方文档 bwa 篇

参考文章：

BWA使用详解：https://www.jianshu.com/p/19f58a07e6f4

https://github.com/lh3/bwa

https://bio-bwa.sourceforge.net/bwa.shtml

https://bio-bwa.sourceforge.net/

BWA介绍

将测序数据比对到参考基因组上，因为精确且高效的读取比对是许多下游分析（如变异检测、RNA-seq分析、结构变异识别等）的基础。

算法：BWA使用Burrows-Wheeler Transform (BWT) 算法来索引参考基因组，这种算法提供了快速且内存效率高的比对。

模式：BWA提供多种比对模式 (BWA-MEM、BWA-SW、BWA-backtrack)，每种模式下有各自的参数可以调整。

数据格式：BWA的输出的数据格式 (SAM/BAM格式) 方便与其它工具使用 (SAMtools、GATK、Picard等) 高度兼容，方便进一步的数据处理和分析。

BWA参数介绍

STEP01: 在进行比对之前需要建索引

bwa index [参考基因组.fasta]

# 实例 
# index/B73_v4_LJS.SNPs.fasta 是对应的 参考基因组 fasta 文件位置
bwa index index/B73_v4_LJS.SNPs.fasta

STEP02: 比对 (Alignment)

三种比对模式 BWA-MEM、BWA-SW、BWA-backtrack

  1. BWA-backtrack（用于较短的读取）
bwa aln [参考基因组.fasta] [reads.fq] > output.sai
bwa samse [参考基因组.fasta] output.sai [reads.fq] > output.sam

# 实例
ref="index/B73_v4_LJS.SNPs.fasta"
bwa aln -t 60 $ref clean/$file.1.fq.gz > $file.1.sai
bwa aln -t 60 $ref clean/$file.2.fq.gz > $file.2.sai

# 使用 bwa sampe 命令将之前用 bwa aln 创建的 .sai 索引文件和原始的 FASTQ 测序文件合并为一个 SAM（Sequence Alignment/Map）格式的输出文件。
# -a 1000 只有当两个配对的 read 在参考序列上距离小于或等于1000碱基对时，它们才会被认为是有效的配对。
bwa sampe -a 1000 -r '@RG\tID:$file\tPL:UNKNOWN\tLB:$file:4\tSM:$file' $ref $file.1.sai $file.2.sai \ clean/$file.1.fq.gz clean/$file.2.fq.gz > bam/$file.sam

需要注意使用BWA-backtrack模式即aln参数时，只会得到sai格式文件，还需要使用bwa samse产出sam文件。

  1. BWA-MEM（用于长读取，且是最常用的模式）
bwa mem -t [线程数] [参考基因组.fasta] [reads1.fq] [reads2.fq] > output.sam

# 实例 -- 改写上面的 bwa aln + bwa sampe 模式
# -R 添加只读组（Read Group）信息到输出的SAM文件
ref="index/B73_v4_LJS.SNPs.fasta"
bwa mem -t 60 -R '@RG\tID:$file\tPL:UNKNOWN\tLB:$file:4\tSM:$file' $ref clean/$file.1.fq.gz clean/$file.2.fq.gz > bam/$file.sam
  1. BWA-SW

BWA-SW 是 BWA 的一个变种，专门用于对长序列进行局部比对。这通常用于对比如 PacBio 或 Nanopore 这样的长读取进行比对

bwa index [参考基因组.fasta]
bwa bwasw [参考基因组.fasta] [长读取.fq] > output.sam
或
bwa bwasw [参考基因组.fasta] [reads1.fq] [reads2.fq] > output.sam

# 实例
# SE - 单端数据
bwa bwasw genome long_read.fq > aln.sam
# PE - 双端数据
bwa bwasw genome read1.fq read2.fq > aln-pe.sam

samtools 处理 SAM/BAM 格式文件

使用 samtools 进行处理，为之后的的 GATK 等工具使用准备

上面的步骤都是为了获得到SAM/BAM文件，之后还需要使用samtools进行处理数据，比如去重复、排序、索引等。这里给出几个使用方法。

samtools view -bS output.sam > output.bam
samtools sort output.bam -o sorted_output.bam
samtools index sorted_output.bam

samtools view -bS bam/$file.sam -o bam/$file.bam -@ 12

  • samtools: 是一个用于处理各种基因组数据格式的工具。
  • view -bS: 将SAM格式转换为更紧凑的BAM格式。
  • -@ 12: 使用12个线程进行操作。
  • -o bam/$file.bam: 输出的BAM文件。

samtools sort --threads 10 -o bam/$file.bam bam/$file.bam

  • sort: 对BAM文件进行排序，这是很多后续分析（例如变异检测）的必要步骤。
  • --threads 10: 使用10个线程进行排序。
  • -o bam/$file.bam: 排序后覆盖原有的BAM文件。

samtools index bam/$file.bam

  • index: 创建BAM文件的索引，这允许更快地随机访问文件，对于一些工具是必需的。

BWA快速使用

这是我在玉米变异定位中具体使用的 BWA 部分的代码

# 因为下面bwa aln并没有用到index文件夹构建的索引，所以可以不用构建索引 -- 直接使用 bwa aln
bwa index index/B73_v4_LJS.SNPs.fasta
ref="index/B73_v4_LJS.SNPs.fasta"

# bwa 比对 alignment
bwa aln -t 60 $ref clean/$file.1.fq.gz > $file.1.sai
bwa aln -t 60 $ref clean/$file.2.fq.gz > $file.2.sai
# bwa sampe 组合 之前 bwa aln 生成的 sai 格式文件和 fq 测序数据 合并为一个 SAM 格式文件
bwa sampe -a 1000 -r '@RG\tID:$file\tPL:UNKNOWN\tLB:$file:4\tSM:$file' $ref $file.1.sai $file.2.sai clean/$file.1.fq.gz clean/$file.2.fq.gz > bam/$file.sam