GWAS走起~（3 bwa）
序列比对的作用：
1、基于相似性的鉴定及功能预测（如果序列保守的话，那么功能也可能会相似）
2、可能在同源进化上有共同祖先（序列相似性是构建演化树的重要依据之一）
与参考序列比对:这方面的工具也很多。(blast、blat、seqmap这些都是老掉牙的了.问什么告诉你,让你知道,我也做公课了,而且,如果有人让你用blast做全基因组,他肯定是不喜欢你,除此,还有一些在线比对软件)
我这里用的是BWA: BWA是用于针对大型参考基因组（例如人类基因组）比对测序读数的程序。它有两个主要组成部分，一个用于读取小于150bp的片段，另一个用于较长的读取片段。
细点说：BWA主要使用三种算法，分别是BWA-backtrack， BWA-SW和BWA-MEM。其实之前短序列比对主要使用的是bwt算法。BWA-backtrack主要用于短点的reads，主要是100bp以下的，而BWA-SW、BWA-MEM主要用于长点的序列，reads长度可以从70bp到1M，如果算法用错就会扩大容错率，导致结果的准确性。
好！开盘。主要分两步盘（1、建立索引；2、比对序列）。
1、 [用参考序列建立索引：bwa index ref.fa
如果参考序列的文件挺大的话，这可能需要等一点时间。
2、序列比对：bwa mem -R '@RG\tID:E1602061\tSM:bar\tLB:library1' -t 4 ref.fna
这是一些参数：了解一下参数
-R 设置reads标头，放到一对引号中，也就是sam文件中的RG部分，为什么要设置RG表头呢，因为同一样品可能包括多个测序结果，来自不同lane，不同文库，或者不同样品的比对结果合并到同一个文件中进行处理，就需要通过RG进行标记区分。RG每个标记用冒号分割键和值，不同标记用 '\t' 分隔。例如'@RG\tID:foo\tSM:bar\tLB:library1'
-M 这个选项也比较重要，因为bwa men在比对的时候，对于一条序列同时比对到基因组不同区域的情况，bwa认为都是最优匹配，但是会与Picard tools不兼容，影响后面GATK检测，这个时候可以设置-M选项，将较短的比对标记为此优，与picard兼容。
-t 设置线程数，多线程可以显著提高比对效率，因为数据比对大，影响还是非常显著的，现在整个行业想进各种方法来提高比对效率，包括提升硬件，在CPU底层优化，改进软件算法等，因为现在一个完整的人基因组比对还需要几个小时，GATK还需要一两天，如果提高计算速度，是非常有意义的。
-T ，给定一个分值，当比对的分值比设置的小时，不输出该比对结果；
-a 将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的reads
-p 智能的pairing，如果不设置，输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired reads进行比对。如果设置-p，则仅以第1个文件作为输入，输入的文件若有2个，则被忽略之
下面在看一下bwasw比对，因为使用了 Smith-Waterman，特别适合长reads比对，因为能够容许更多的gap情况。帮助文档中默认提示的就是query.fa，一般是长序列。帮助中提示对于长的Illumina 数据, 454，Sanger拼接的contigs或者 fosmids，BACs等，使用默认参数就行，如果是2010年之前Pacbio，设置'-b5 -q2 -r1 -z10'。
bwa bwasw genome long_read.fq > all.sam
bwa bwasw genome read1.fq read2.fq > all.sam
反正不同的的情况,要用不同的算法,自己要看一下官网的介绍.
这步完成会有一个sam文件，这个文件装载着咱们的比对信息，也是bwa软件最终给出的包含比对信息的文件,后面还要将这些比对信息更加凸显化。哦,等一等,因为sam文件是多维数据的文本文件,所以文件很大,处理起来也就比较费事,但是前辈们就将它转化为小点的文件(bam文件)进行处理,就方便多了,下面我们就学习一下这个文件转换工具samtools吧。
