BBQ（生信基础问题16）-BWA算法原理及软件实用

  从这个BBQ开始，我们进阶到实战分析，从转录组开始，一点点的探寻原理，一步步实现！！！！！
--------------------------------------------------分割线----------------------------------------------

我们来思考以下几个主要的大问题：

1：我们利用RNA-seq这个手段可以解决什么问题？
2：RNA-seq的基本原理是什么？
3：我们如何去看这个RNA-seq出的结果？
4：在这个过程中我们运用到的软件有哪些？什么原理？

这些问题都比较的大，我们可以先Mark一下，然后我们在学习过程中，一点点的去解决！

--------------------------------------------分割线------------------------------------------------------
之前一直在说序列比对的问题，以及一些数据下机的质控等等。

一：我们来想一想，入门拿到了这些数据数据之后，我们如何去处理这些数据？

  1.1:根据我们之前的学习的知识，肯定是需要对下机的数据进行质量评估，那么用到的软件是什么？

数据下机后的常规操作

  1.2：拿到了这些clean data我们应该怎么办？序列比对！那么和什么序列比对，比对的目的是什么？用什么软件比对？

  二代测序技术或者说是高通量测序技术产生的数据有几个明显的特征：1是测序数据量大，2是序列比较短，3是测序错误一般可以控制在0.1%左右。那么如何快速将这些测序数据比对到参考基因组是一个比较重要的问题。

  对于数据分析来说，输入常常是10^7 甚至更多的reads，而且是要在全基因组上寻找定位，比如人的基因组有3Gbp大！所以如果不优化算法，估计mapping这个问题就要等到地老天荒。

  将reads map到参考基因组，现在主流的序列比对的方法有：bowtie, bowtie2 or bwa。Bowtie2在Bowtie的基础上进行了一次更新，但是核心原理基本没有改变，只是优化了下游部分的比对策略，在bowtie的基础上允许了deletion的出现。但是他们的核心都是基于BWT算法。

1.3:什么是BWT算法？

  BWT (Burrows–Wheeler_transform)数据转换算法：个人理解就是对数据压缩建立索引，然后方便进行后续的比对。如何压缩，如何转换？
我们以S='acaacg'这个字符串为例子

  1.3.1在原本的字符串后面添加一个额外的字符：美刀符号.

那么原来的字符就变成了acaacg$，
注意：新加入的这个字符一定要比原序列中的字符排序上“小”。

  1.3.2 把这个新生成的字符串进行矩阵变化,说白了就是把美刀符号往前移动了一位，这样子就得到了一个矩阵。如下图

新生成的矩阵

  1.3.3对矩阵进行排序：根据字符的“大小”进行排序，还记得谁是第一个吗？（必须是那个含有美刀符号的序列呀~）

排序后的序列，注意看第一排！

  1.3.4：取到最后一列和第一列，分别为F列（First）和L列（Last）.

F列和L列

  1.3.5进行回溯啦。我们需要找到L列的最小的那个字符（也就是$符号）。找到它在F列的位置。

回溯过程1
我们找到了它在L列的位置，它对应的是哪个字符呢（g），所以我们确定了第一个顺序g.
回溯2
这个g对应了F连的哪个位置呢？看箭头~
回溯3
这个g对应了哪个位置呢？（我们可以轻松的找到是C）
回溯4
但是问题来了，原序列有2个c怎么办？这应该找哪个对应呢？
按照顺序来：g回溯到的是第二个c，那么也得对应的是第二个c
回溯5（图丑了点，凑合看吧）
按照这个方法，我们可以得到'$gcaaca'，于是S='acaacg'

  1.3.6需要注意的问题：

1）首先，F和L中同行的两个字符在原字符串中必然是相邻的，并且L中的字符必然在F中的字符的前面。

2）在L中有时会有同一个字符的多次出现，它们在L中的顺序和F中的顺序是相同的。比如，在L中的第一个出现的a，它是T中的第三个a在F中的第一个出现的a也是T中的第3个a。

根据这个性质，那么我们可以得到一个重要的推论，就是在L中的每一个字符都可以唯一的在F中找到，且在T中的位置是与之相对应的，即同一个字符的在T中的位置仍然保持不变，FM-index就是就是利用此规律进行查询的.

BBQ问题：

1. 为什么FASTQ文件的快速比对需要建立index

因为如果直接进行序列比对的话，效率慢耗时长。建立index可以更快速的进行比对。

2. 如果我从1个网站上下载的是1个物种的参考转录组的序列，其中包含了A,U,C,G碱基，我的FASTQ为该物种转录组测序的结果，用A,G,T,C，4种碱基来表示。那么需不需要在建立index之前把参考转录组中的U全部都换成T？>>肯定要进行转换，但是我觉得这个问题略白，因为进行RNA测序的时候，都是进行了转换变成了cDNA，所以无论是参考基因组还是下机的数据都是ATCG，但是如果有U的话肯定要转换成T。这只能自己写脚本了

二：软件实战

BBQ3:请在Linux环境下，下载human genome 19参考基因组的1号染色体序列；并使用bowtie 建立index。

1:wget ./test http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr1.fa.gz
2:gzip -d ./test/chr1.fa.gz
3:bowtie2-build ./test/chr1.fa ./test/chr_bwt2_index 

Reference:
1:生物信息学100个基础问题 —— 第16题 高通量测序的回贴问题 I
2:踏踏实实做技术：BWA，Bowtie，Bowtie2的比对算法推导
3: DNA比对算法：BWT
4:BWA的设计思想2