WES2Neoantigen Pipeline

Part 2 BWA

一、建立索引

$ bwa index -a bwtsw human.fa

-a [is|bwtsw] : 输入构建Index的算法。is算法快速简单，是默认选项，但是不能用于基因组大于2GB的数据库。bwtsw适用于大基因组。

-p STR：输出的数据库的前缀。默认与输入的文件名一致。

完成之后，你会看到类似如下几个以human.fasta为前缀的文件：

.

├── human.fasta.amb

├── human.fasta.ann

├── human.fasta.bwt

├── human.fasta.pac

└── human.fasta.sa

这些就是在比对时真正需要被用到的文件(因此很多网站在提供参考序列的时候，连索引文件也一起提供，下载索引比自己建立索引文件要快很多)。这一步完成之后，我们就可以将read比对至参考基因组了。

二、开始比对

BWA mem比对模块

该算法先使用MEM（maximal exact matches)进行seeding alignments,再使用SW算法（affine-gap Smith-Waterman)算法进行seed的延伸。BWA-MEM算法执行局部比对和剪接性。可能会出现query序列的多个不同的部位出现各自的最优匹配，导致reads有多个最佳匹配位点。

Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]

bwa mem -t 4 -M genome read1.fq read2.fq > aln-pe.sam

-t 使用的线程数

-p 若无此参数，输入文件只有1个，则进行单端比对，输出文件有两个，则作为paired reads进行比对。若加入此参数，则仅以第一个文件作为输入，该文件必须是read1.fq和read2.fq进行reads交叉的数据

-M 加入此参数用于将shorted split hits标记为次优，有利于兼容Picard.

-k 最小种子匹配长度，影响比对速度，默认为19

-w gap最大长度，超过后将不会被匹配，默认为100

-d 影响停止延伸的一个参数。(Off-diagonal X-dropoff (Z-dropoff). Stop extension when the difference between the best and the current extension score is above |i-j|*A+INT, where i and j are the current positions of the query and reference, respectively, and A is the matching score. Z-dropoff is similar to BLAST’s X-dropoff except that it doesn’t penalize gaps in one of the sequences in the alignment. Z-dropoff not only avoids unnecessary extension, but also reduces poor alignments inside a long good alignment. [100])

-r 一个关键的启发式搜索参数，较大的数值意味着较少的种子，比对速度也更快，默认1.5，即1.5倍的最小种子长度

-c INT  Discard a MEM if it has more than INT occurence in the genome. This is an insensitive parameter. [10000]

-A INT  匹配得分，默认1

-B INT 错配惩罚，默认4. The sequence error rate is approximately: {.75 * exp[-log(4) * B/A]}. [4]

-O INT  Gap open penalty. [6]

-E INT  Gap extension penalty. A gap of length k costs O + k*E (i.e. -O is for opening a zero-length gap). [1]

-L INT  Clipping penalty. When performing SW extension, BWA-MEM keeps track of the best score reaching the end of query. If this score is larger than the best SW score minus the clipping penalty, clipping will not be applied. Note that in this case, the SAM AS tag reports the best SW score; clipping penalty is not deducted. [5]

-U INT  Penalty for an unpaired read pair. BWA-MEM scores an unpaired read pair as scoreRead1+scoreRead2-INT and scores a paired as scoreRead1+scoreRead2-insertPenalty. It compares these two scores to determine whether we should force pairing. [9]

-p Assume the first input query file is interleaved paired-end FASTA/Q. See the command description for details.

-R 设置reads标头，放到一对引号中，也就是sam文件中的RG部分。RG每个标记用冒号分割键和值，不同标记用 '\t' 分隔。例如'@RG\tID:foo\tSM:bar\tLB:library1' (为什么要设置RG表头呢，因为同一样品可能包括多个测序结果，来自不同lane，不同文库，或者不同样品的比对结果合并到同一个文件中进行处理，就需要通过RG进行标记区分)。1) ID，这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；2) PL，指的是所用的测序平台，这个信息不能随便写！，特别是当我们需要使用GATK进行后续分析的时候，更是如此！在GATK中，PL只允许被设置为：ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这几个信息。3) SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本；4) LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB；

-T 给定一个分值，当比对的分值比设置的小时，不输出该比对结果。

-a 将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的 reads。

-C  Append append FASTA/Q comment to SAM output. This option can be used to transfer read meta information (e.g. barcode) to the SAM output. Note that the FASTA/Q comment (the string after a space in the header line) must conform the SAM spec (e.g. BC:Z:CGTAC). Malformated comments lead to incorrect SAM output.

-H  Use hard clipping ’H’ in the SAM output. This option may dramatically reduce the redundancy of output when mapping long contig or BAC sequences.

-M 这个选项也比较重要，因为bwa men在比对的时候，对于一条序列同时比对到基因组不同区域的情况，bwa认为都是最优匹配，但是会与Picard tools不兼容，影响后面GATK检测，这个时候可以设置-M选项，将较短的比对标记为次优，与picard兼容。

-v INT  Control the verbose level of the output. This option has not been fully supported throughout BWA. Ideally, a value 0 for disabling all the output to stderr; 1 for outputting errors only; 2 for warnings and errors; 3 for all normal messages; 4 or higher for debugging. When this option takes value 4, the output is not SAM. [3]

Eg：

bwa index ref.fna

bwa mem -R '@RG\tID:E1602061\tSM:bar\tLB:library1' -t 4 ref.fna reads.1.fq.gz reads.2.fq.gz  >result.sam

注意：

1、安装时需要依赖zlib；

2、建库时注意选择正确的算法；

3、最新版本bwa，需要使用-R选项。

测试实例：

/pub5/xiaoyun/BioSoftware/BWA/bwa-0.7.12/bwa mem -R '@RG\tID: ERR863675\tPL:ILLUMINA\tSM:ERR863675\tLB:library1' -t 20 /IJob/J33/HLA_example_temp/temp/Sequence/reference.hg19.fa /pub6/temp/WLY/Trimmomatic_Result/out.ERR863675_1.fastq /pub6/temp/WLY/Trimmomatic_Result/out.ERR863675_2.fastq > /pub6/temp/WLY/BWA_Result/AfterQC.sam

system(paste(BWA,"mem","-t 20","-R '@RG\tID:Test\tSM:Test\tPL:ILLUMINA\tLB:Library'",Reference,"/IJob/J33/Data/SRP102119_BC_12P22M34NSamples_WES/fastq/_AfterFASTQC_/AlignmentTest/SRR5357679_1P.fastq.gz","/IJob/J33/Data/SRP102119_BC_12P22M34NSamples_WES/fastq/_AfterFASTQC_/AlignmentTest/SRR5357679_2P.fastq.gz","1>/IJob/J33/Data/SRP102119_BC_12P22M34NSamples_WES/fastq/_AfterFASTQC_/AlignmentTest/SRR5357679.sam","2>/IJob/J33/Data/SRP102119_BC_12P22M34NSamples_WES/fastq/_AfterFASTQC_/AlignmentTest/SRR5357679.align.log",sep=" "))

Ref:

1、知乎

2、知乎

3、CSDN

4、简书