NGS数据比对之BWA

BWA主要是将reads比对到大型基因组上，主要功能是：序列比对。首先为大型参考基因组建立索引，然后将reads比对到基因组。特点是快速、准确、省内存。由三种类似算法组成：BWA-backtrack，BWA-SW和BWA-MEM。首推BWA-MEM。

三种算法的适用范围

BWA-backtrack：reads长度<70bp时，推荐本算法，建议输入reads长度 < 100bp。
BWA-SW：在reads具有频繁的gap时，比对更敏感，推荐本算法。reads长度一般为70bp-1Mbp，支持long-reads，split alignment。
BWA-MEM(首推)：在reads长度在70bp-1Mbp范围时，推荐本算法(除了上面两种情况)。支持long-reads，split alignment。

BWA使用说明 使用手册

语法
bwa index ref.fa #首先建立基因组索引
bwa mem ref.fa reads.fq > aln-se.sam # 调用BWA-MEM
bwa mem ref.fa read1.fq read2.fq > aln-pe.sam # 调用BWA-MEM
bwa aln ref.fa short_read.fq > aln_sa.sai # 调用BWA-backtrack
bwa samse ref.fa aln_sa.sai short_read.fq > aln-se.sam # 调用BWA-backtrack
bwa sampe ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam # 调用BWA-backtrack
bwa bwasw ref.fa long_read.fq > aln.sam # 调用BWA-SW
注意：BWA输入的是fastq/fq的原始测序数据。

bwa的使用需要两个 输入文件 ：
Reference genome data 和 Short reads data

Step1.建立索引

根据reference genome data 建立 Index File

bwa index -a bwtsw reference.fa   

##bwa index [-p prefix] [-a algoType] <in.db.fasta>
bwa index -a is TP53.fna

构建索引时需要注意的问题：bwa构建索引有三种算法，三种算法都是基于BWT的，这三种算法通过参数
-a is 、－a div和-a bwtsw进行选择。其中-a bwtsw对于短的参考序列是不工作的，必须要大于等于10Mb；-a is（效果和-a div是一样的）是默认参数，这个参数不适用于大的参考序列，必须要小于等于2G。

Step2.比对 bwa mem

image.png

module add bwa
module add broad_hg38/bwa
SAMPLE=S001P1
RESULT_ROOT="sunyq"

echo "bwa starts at `date`"
bwa mem -R "@RG\\tID:${SAMPLE}\\tCN:NovoGene\\tLB:SureSelectHumanAllExonV6r2\\tPL:Illumina\\tSM:${SAMPLE}" \
        -K 100000000 -v 3 -t 2 -Y \
          $GENOME_LOCAL/$GENOME_NAME \
          S001P1_1.fq.gz \
          S001P1_2.fq.gz \
       -o ${RESULT_ROOT}/${SAMPLE}_bwa_hg38.sam

echo "bwa ends at `date`"

需要参考基因组以及原始fastq文件，输出sam文件格式

sam文件解读：

read包含:

image.png

FLAG字段中的定义为:

image.png

BWA生成以下可选字段。以“X”开头的标签是特定于BWA的。

image.png

使用bwa完成比对后，用samtools完成BAM格式转换、排序并标记PCR重复序列。

例子：利用 SRR1770413 ，见WES数据处理之寻找突变GATK

bwa index -a is  E.coli_K12_MG1655.fa  #创建参考基因组索引

module add bwa

bwa mem -R "@RG\tID:foo\tPL:illumina\tSM:E.coli_K12" -t 2 /home/linjc/jianshu/WES/E.coli_K12_MG1655.fa \
          /home/linjc/jianshu/WES/SRR1770413/E_Coli_CGATGT_L001_R1.fastq.gz \
          /home/linjc/jianshu/WES/SRR1770413/E_Coli_CGATGT_L001_R2.fastq.gz \
       -o /home/linjc/jianshu/WES/E_Coli_CGATGT_bwa.sam

samtools view -Sb E_Coli_CGATGT_bwa.sam > E_Coli_CGATGT_bwa.bam  #转换为bam文件

#排序
samtools sort -@ 4 -m 4G -O bam -o /home/linjc/jianshu/WES/E_coli_K12.sorted.bam /home/linjc/jianshu/WES/E_Coli_CGATGT_bwa.bam   

#标记重复PCR  module add gatk/4.2.6.1
gatk MarkDuplicates -I /home/linjc/jianshu/WES/E_coli_K12.sorted.bam -O /home/linjc/jianshu/WES/E_coli_K12.sorted.markdup.bam -M /home/linjc/jianshu/WES/E_coli_K12.sorted.markdup_metrics.txt

#创建比对索引文件
samtools index /home/linjc/jianshu/WES/E_coli_K12.sorted.markdup.bam

image.png

以上为BWA比对流程，排序，标记等