生信基础知识

NGS数据比对之BWA

2022-05-18  本文已影响0人  忙碌的dog

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比对流程,排序,标记等

上一篇下一篇

猜你喜欢

热点阅读