使用bwa进行序列比对

 

001、

bwa mem -t 4 -k 32 -M -R "@RG\tID:name\tSM:name\tPL:illumina\tLB:name\tPU:name" reference.fna sm.clean.1.fastq.gz sm.clean.2.fastq.gz | samtools view -Sb - > sm.bam 

 

mem：mem比对算法

-t：指定线程数

-k： （这个参数可以不设置）最小的种子长度（minimum seed length），bwa采用seed-and-extend的策略，在seed阶段，bwa取read的碱基片段在reference上进行精确匹配，并选择满足一定长度要求和匹配次数的read片段作为seed，这个阶段算法的核心是基于FM-index的精确匹配； 在extend阶段，bwa利用smith-waterman算法将seed在read和reference上向两边延伸比对（容忍gap），进而找到整个read在reference上符合条件的全局匹配。

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

-R：Read Group的字符串信息，这是一个非常重要的信息。 

(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这几个信息。基本上就是目前市场上存在着的测序平台，当然，如果实在不知道，那么必须设置为UNKNOWN，名字方面不区分大小写。

(3) SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本。

 

(4) LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB；

除了以上这四个之外，还可以自定义添加其他的信息，不过如无特殊的需要，对于序列比对而言，这4个就足够了。这些信息设置好之后，在RG字符串中要用制表符（\t）将它们分开。

 

管道后接samtools：

-S：表示读入文件为sam格式

-b：表示输出文件为bam格式

 。

 

参考：https://developer.aliyun.com/article/1242282

。