Ⅰ BWT比對演算法
BWT演算法在多款序列比對軟體（BWA，bowtie）中都有涉及，那麼對於RNA-seq的2代數據，一般建庫長度是單端300bp，雙端各150bp左右。
對於兩個序列進行比對，即pairwise alignment，我們可以按比對方式分為全局比對（NW演算法）和局部比對（SW演算法）：
比方說SOAP：
它把參考基因組分成一個個小片段，並且將這些小片段建立成hash關系，也就是說我在序列比對的時候，只要我獲取到reads的內容，我就可以match（翻譯為完全比對？）到參考基因組上，並獲得其位置信息。
這里介紹的是BWA，bowtie。而這兩款軟體除了要對基因組建立index以外，還採用的是BWT演算法來實現比對。
step2 ：
之後，「$」符號向後平移一個單位，得到 :
「$」符號再向後平移一個單位:得到:
以此類推，最終得到：
這樣的序列的排序表，接下來按照字母順序進行排序（$,A,C,G,T）有：
這個矩陣我們稱之為轉換矩陣
我們根據F列和L列元素的相對位置就可以推導出該序列的全部信息：
首先，我們單獨取出F列和L列
在F列和L列中找到「$」，連接它們，那麼在水平位置對應於的G（藍色箭頭）
這個G就是「$」前面的元素:
接下來在F列中找到元素G：
此時按照上述方法，才F列中找到元素C，由於L列中的元素C為第二個C（上面G後面的C視為第一個C），所以應該對於F列的第二個C
接下來L列第三個A對應F列第三個A，因此以此類推就可以最終還原（最終對應於L列的$)：
比方說我現在有一條短序列CAA要比對到ACAACG上，首先說明的是在比對過程中，CAA是反向進行比對的，即A->A->C
由於CAA第一個比對的元素是A，在F列中有三個A，因此我們分類討論
選擇第一個A：
此時比對結果為ACA，不符合
選擇第二個A：
結果為$，不符合
選擇第三個A：
此時比對結果為CAA，符合原序列
在現實中，CAA可以想像為fq文件reads，ACAACG可以想像為參考基因組，並且在比對中會存在mismatch和gap這種情況的發生，因此軟體會對每一種比對的情況進行打分，擇優處理
目前 bowtie 是沒有gap這一項設定的
參考： https://www.bilibili.com/video/av15743137/
劉小樂哈佛大學課程： https://www.bilibili.com/video/BV1De411s71p
https://www.bilibili.com/video/BV1d4411E7uS?from=search&seid=7482313674070699968
Ⅱ 10.12 bwa使用 安裝文件路徑與使用 sh許可權
我們這里將用於流程構建的BWA就是其中最優秀的一個，它將BW(Burrows-Wheeler)壓縮演算法和後綴樹相結合，能夠讓我們以較小的時間和空間代價，獲得准確的序列比對結果。
別人的已安裝文件打包傳遞後使用：
1、連接伺服器
2、家目錄下，Users，ls -all (或者 打開/etc里profile）
3、vim .bashrc
```
export PATH="your path:$PATH"
```
4、source .bashrc
5、在全局中輸入bwa測試是否可以使用
自己安裝的軟體路徑都已經配置好了，所以不需要這樣操作就可以直接使用。
BWA MEM比對模塊是有一定適用范圍的：它是專門為長read比對設計的，目的是為了解決，第三代測序技術這種能夠產生長達幾十kb甚至幾Mbp的read情況。一般只有當read長度≥70bp的時候，才推薦使用，如果比這個要小，建議使用BWA ALN模塊。
BWA ALN模塊
1: 建立 Index
根據reference genome data(e.g. reference.fa) 建立 Index File
bwa index -a bwtsw reference.fa
(bwa index 可以查看更多使用方法）
step 2: 尋找 SA coordinates
如果是pair-end 數據（leftRead.fastq和rightRead.fastq）兩個文件分別處理
bwa aln reference.fa leftRead.fastq > leftRead.sai
bwa aln reference.fa rightRead.fastq > rightRead.sai
bwa aln reference.fa singleRead.fastq > singleRead.sai
如果希望多線程運行，在其中加入 -t這個參數，另外-f這個參數可以指定結果輸出文件，如:
bwa aln -c -t 3 -f leftreads.sai reference.fa leftreads.fastq
step 3:轉換SA coordinates輸出為sam
如果是pair-end數據
bwa sampe -f pair-end.sam reference.fa leftRead.sai rightRead.sai leftRead.fastq rightread.fastq
如果是single reads數據
bwa samse -f single.sam reference.fa single.sai single.fastq
https://www.jianshu.com/p/859c0345624c
bwa mem [options] <idxbase> <in1.fq> [in2.fq]
bwa mem -t 4 -R'@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name'/path/to/human.fasta read_1.fq.gz read_2.fq.gz > sample_name.sam
-t，線程數，我們在這里使用4個線程；-R 接的是 Read Group的字元串信息，這是一個非常重要的信息 ，以@RG開頭，它是用來將比對的read進行分組的。不同的組之間測序過程被認為是相互獨立的，這個信息對於我們後續對比對數據進行錯誤率分析和Mark plicate時非常重要。在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，名字方面不區分大小寫。如果你在分析的時候這里沒設置正確，那麼在後續使用GATK過程中可能會碰到類似如下的錯誤。
(3) SM，樣本ID，同樣非常重要，有時候我們測序的數據比較多的時候，那麼可能會分成多個不同的lane分布測出來，這個時候SM名字就是可以用於區分這些樣本。
(4) LB，測序文庫的名字，這個重要性稍微低一些，主要也是為了協助區分不同的group而存在。文庫名字一般可以在下機的fq文件名中找到，如果上面的lane ID足夠用於區分的話，也可以不用設置LB；
除了以上這四個之外，還可以自定義添加其他的信息，不過如無特殊的需要，對於序列比對而言，這4個就足夠了。這些信息設置好之後， 在RG字元串中要用製表符（\t）將它們分開 。
sam轉換bam（sam特殊二進制）
bwa mem -t 4 -R'@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name'/path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - > sample_name.bam
time$bwamem -t 40 -M -Y -R"@RG\tID:$RGID\tPL:ILLUMINA\tSM:$sample"$reference\$outdir/cleanfq/${fq_file_name}.paired.1.fq.gz$outdir/cleanfq/${fq_file_name}.paired.2.fq.gz |$samtoolsview -Sb - >$outdir/bwa/${sample}.bam && \echo"** BWA MEM done **"&& \
-M 將 shorter split hits 標記為次優，以兼容 Picard』s markDuplicates 軟體。
-Y Use soft clipping CIGAR operation for supplementary alignments. By default, BWA-MEM uses soft clipping for the primary alignment and hard clipping for supplementary alignments.
鏈接：https://www.jianshu.com/p/d6cc153be28d
設置許可權 chmod 777 ./*.sh
運行
以上，我們就完成了read比對的步驟。接下來是排序：
-@，用於設定排序時的線程數；-m，每個線程排序時最大的內存消耗；-O 指定輸出為bam格式；-o 是輸出文件的名字.
在排序完成之後我們就可以開始執行去除重復（准確來說是 去除PCR重復序列 ）的步驟了。
把參數REMOVE_DUPLICATES設置為ture，那麼重復序列就被刪除掉，不會在結果文件中留存。我比較建議使用第一種做法，只是標記出來，並留存這些序列，以便在你需要的時候還可以對其做分析。
這一步完成之後，我們需要為sample_name.sorted.markp.bam創建索引文件，它的作用能夠讓我們可以隨機訪問這個文件中的任意位置，而且後面的「局部重比對」步驟也要求這個BAM文件一定要有索引，命令如下：
$ samtools index sample_name.sorted.markp.bam
接下來是局部區域重比對，通常也叫Indel局部區域重比對。有時在進行這一步驟之前還有一個merge的操作，將同個樣本的所有比對結果合並成唯一一個大的BAM文件【注】，merge的例子如下：
$ samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]
【注意】之所以會有這種情況，是因為有些樣本測得非常深，其測序結果需要經過多次測序（或者分布在多個不同的測序lane中）才全部獲得，這個時候我們一般會先分別進行比對並去除重復序列後再使用samtools進行合並。