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
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包含:
FLAG字段中的定义为:
BWA生成以下可选字段。以“X”开头的标签是特定于BWA的。
使用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
以上为BWA比对流程,排序,标记等