一、轉錄組還是基因組？

map常用的工具有bowtie/bowtie2, BWA,SOAP1/SOAP2等。這個問題又會被分成兩個問題，是基因組測序（DNA-seq）還是轉錄組測序(mRNA-seq)。其中的區別是對於真核生物而言，mRNA序列與DNA序列並不完全相同，在經歷了後剪下之後，成熟的mRNA可能是原基因的一部分，甚至順序及個別鹼基會產生變化。如果是mRNA測序，那map工作就會在DNA測序map的基礎上再多一步，map到轉錄組上去。所以最為流行的做法是，（使用BWA來進行ChIP-seq測序）使用bowtie來map DNA測序，使用tophat來map RNA測序。實際上，tophat是通過呼叫bowtie來完成工作的。而tophat1和tophat2的差別最主要的就是呼叫了bowtie1還是bowtie2。當然如果你只安裝了bowtie1的話，tophat2也可以呼叫它

二、bowtie有1和2的差別：

1，bowtie1出現的早，所以對於測序長度在50bp以下的序列效果不錯，而bowtie2主要針對的是長度在50bp以上的測序的。
2，Bowtie 2支援有空位的比對
3，Bowtie 2支援區域性比對，也可以全域性比對
4，Bowtie 2對最長序列沒有要求，但是Bowtie 1最長不能超過1000bp。
5. Bowtie 2 allows alignments to [overlap ambiguous characters] (e.g. `N`s) in the reference. Bowtie 1 does not.
6，Bowtie 2不能比對colorspace reads.
…………

三、選不選bowtie：

Bowtie 2 在比對大的基因組的時候，工作效果更好一些，如果你是比對兩個特別大的基因組，可以考慮MUMmer，如果你比對的是兩個特別相近且比較小的基因組，可以用[NUCmer], [BLAT], or [BLAST]. 但是Bowtie 2不能比對colorspace reads.

四、Bowtie的簡介：

Bowtie是一個超級快速的，較為節省記憶體的短序列拼接至模板基因組的工具。它在拼接35鹼基長度的序列時，可以達到每小時2.5億次的拼接速度。
Bowtie並不是一個簡單的拼接工具，它不同於Blast等。它適合的工作是將小序列比對至大基因組上去。它最長能讀取1024個鹼基的片段。模板最小尺寸不能小於1024鹼基，而短序列最長而不能超過1024鹼基。換言之，bowtie非常適合下一代測序技術。
在 使用bowtie前，需要使用bowtie-build來構建比對模板。

Bowtie設計思路是：1）短序列在基因組上至少有一處最適匹配， 2）大部分的短序列的質量是比較高，3）短序列在基因組上最適匹配的位置最好只有一處。這些標準基本上和RNA-seq, ChIP-seq以及其它一些正在興起的測序技術或者再測序技術的要求一致。

使用tophat的基本步驟是：
理解bowtie
使用tophat來mapping reads。其命令常見的形式為：
/path-to-bowtie-programs/tophat -o -p cpu /path-to-genome-index/prefix fastq file 
For example:
/programs/tophat -o TophatOutput/ -p 8 /programs/indexes/hg19 Experiment1.fastq
Paired-end Example:
/programs/tophat -o TophatOutputPE/ -p 8 /programs/indexes/hg19 Experiment1.r1.fastq Experiment1.r2.fastq
可以發現，其很多引數是同bowtie是一樣的。但是它有幾個重要引數需要了解：
--library-type
 用於生成RNA-seq的library。最常見的是使用fr-unstranded，兩條鏈都考慮。
-G   用於加註transcriptome資訊。

五、Bowtie 2對example資料夾中Lambda的處理

1,對lambda_virus.fa建立索引值

cd /sam/bowtie2/example/myindex
 /sam/bowtie2/bowtie2-build /sam/bowtie2/example/reference/lambda_virus.fa lambda_virus
#將會生成以lambda_virus開始名字，字尾為`.1.bt2`, `.2.bt2`, `.3.bt2`, `.4.bt2`,
`.rev.1.bt2`, and `.rev.2.bt2`. 的六個檔案；bowtie2-build 可以處理來自[UCSC], [NCBI],[Ensembl]的fasta檔案，當處理多個fasta檔案的時候，可以用逗號分開處理

2.1 針對unpaired reads的比對

/sam/bowtie2/bowtie2 -p cpu -x lambda_virus -U /sam/bowtie2/example/reads/reads_1.fq -S eg1.sam

2.2針對Paired-end reads的比對

（檔案的格式通常為e.g. `-1 flyA_1.fq,flyB_1.fq`）
/sam/bowtie2/bowtie2 -p cpu -x lambda_virus -1 /sam/bowtie2/example/reads/reads_1.fq -2 /sam/bowtie2/example/reads/reads_2.fq -S eg2.sam
-p跟cpu的個數；
-x跟的是建立的索引，不用跟名字後面的字尾；
-1，-2分別表示Paired-end reads；
輸出格式預設為SAM。
可以指定為多種其它格式，比如說BAM（SAM的二進位制檔案）等等

問題：

這裡的reads用的是clean reads，還是raw reads,我覺得應該是clean reads,如果是clean reads 的話，就得先對raw reads進行trim,但是trimmomatic進行trim後生成的可是4個檔案（因為Illumina測序出來paired end 就有正反兩個raw reads，trim後就包括了能配對的兩個正反reads,和不能配對的兩個正反reads），這樣的話，我用哪個呢，就用僅僅能配對上的兩個clean reads?,如果這樣的話，我的拼接過程中可是都用了的，要不分開去匹配，但是分開匹配的，最後的結果又怎麼統計呢？
先用bowtie處理paired的，然後再用-u處理unpaired 得到兩個sam,再cat ,然後再用samtools sort後來統計結果。
我誤以為這個bowtie既可以處理pair的，也可以同時處理unpair的，結果出現上面的問題。實際上應該是bowtie2處理的 paired reads 在的兩個檔案，必須有同樣多的 reads數，並且都是一一對應的；如果一對不能都比上 會找其中一條來比對上；一個如果不能比對到基因組，那麼使用 –no-mixed 引數，則不會對與之匹配的reads進行比對了。總之核心思想就是必須處理一對，如果一對中有一條未能匹配上，你可以選擇是否保留這一對的結果，–no-mixed就是你的選擇。

2.3區域性比對的例子

用[local alignment] 來比對更長的reads included
 $BT2_HOME/bowtie2 --local -x lambda_virus -U $BT2_HOME/example/reads/longreads.fq -S eg3.sam

問題：這裡的的local跟之前的有啥區別呢？
可以加入-p 12來調動我的12個cpu,速度會更快
檢視生成的eg1.sam檔案
head eg1.sam

2.4其他引數的說明：

-U 後面接的是unpaired reads
 -S 輸出檔案，預設的是stdout，在終端顯示，也可以搞個檔案來裝輸出的結果
輸入檔案選項
 -q Reads 為FASTQ格式（`.fq` or `.fastq`.），此為預設的格式
 --qseq Reads 為QSEQ檔案（end in `_qseq.txt`）. 
 -f Reads 為FASTA檔案（`.fa`, `.fasta`, `.mfa`, `.fna` or similar），如果設定了`-f`,結果相當於`--ignore-quals`也被設定，就是不衡量質量了
 -r Reads為每行就是一條序列，沒有read的名字，沒有質量，出來的結果也相當於`--ignore-quals被設定了
 -c read 序列通過命令列輸入，包含`--ignore-quals`.
 -s/--skip 跳過前面的多少個鹼基再來比對
 -u/--qupto 僅僅比對前面的幾個鹼基，然後停止比對
 -5/--trim5 在比對之前去掉5’端的幾個鹼基（預設的為0）
 -3/--trim3 同理同上
 --phred33 輸入檔案的質量為33
 --phred64 同理同上，在我之前的部落格中有記錄http://blog.sina.com.cn/s/blog_670445240101kq45.html

Bowtie提供了兩個引數（–phred33和–phred64）來指定計分系統，預設是前者。而且，bowtie會根據輸入的具體資料來識別輸入資料實際採用哪套計分系統，使用者不必指定。
針對Phred+33計分系統
!-% 33-37 罰分為0。
&-/ 38-47 罰分為10。
0-9 48-57 罰分為20。
:-I-~ 58-73-126罰分為30。Phred+33實際分數最高到73
針對Phred+64計分系統
!-% 64-68 罰分為0。
&-/ 69-78 罰分為10。
0-9 79-88 罰分為20。
:-I-~ 89-104-126罰分為30。Phred+64實際分數最高到104

--solexa-quals 預設的是關閉的，不解釋
 --int-quals
`--end-to-end` 模式所涉及的幾個引數
 --very-fast Same as: `-D 5 -R 1 -N 0 -L 22 -i S,0,2.50`
 --fast Same as: `-D 10 -R 2 -N 0 -L 22 -i S,0,2.50`
 --sensitive Same as: `-D 15 -R 2 -L 22 -i S,1,1.15` (default in `--end-to-end` mode)
 --very-sensitive Same as: `-D 20 -R 3 -N 0 -L 20 -i S,1,0.50`
`--local` 模式下的幾個引數
 --very-fast-local Same as: `-D 5 -R 1 -N 0 -L 25 -i S,1,2.00`
 --fast-local Same as: `-D 10 -R 2 -N 0 -L 22 -i S,1,1.75`
 --sensitive-local Same as: `-D 15 -R 2 -N 0 -L 20 -i S,1,0.75` (default in `--local` mode)
 --very-sensitive-local Same as: `-D 20 -R 3 -N 0 -L 20 -i S,1,0.50`

比對的結果引數：

-N 在種子比對的過程中容許的錯配數，一般可以設定為0或者1，設的越高，越慢，但是可以增加比對的靈敏性，預設的為0
-L 設定的種子字串的長度，越小，越靈敏，但是越慢，預設的引數設定參考上面的。
-i 怎麼設定種子，詳見官網說明書
……
其他的引數具體看官網說明書

六、使用SAMtools/BCFtools來接下來分析

SAMtoolsis是分析SAM和BAM檔案的工具
BCFtools分析突變和操作VCF和BCF檔案的工具
例子：
1，先用bowtie2生成sam檔案
$BT2_HOME/bowtie2 -x $BT2_HOME/example/index/lambda_virus -1 $BT2_HOME/example/reads/reads_1.fq -2 $BT2_HOME/example/reads/reads_2.fq -S eg2.sam
-x 表示建立的索引的字首名字 -1 -2 後面分別的是雙末端測序的reads –S為輸出的結果
2，samtools 將SAM檔案轉化為BAM檔案
samtools view -bS eg2.sam > eg2.bam
3，用samtools sort將BAM檔案進行排序。
samtools sort eg2.bam eg2.sorted
4，尋找突變通過VCF格式
samtools mpileup -uf $BT2_HOME/example/reference/lambda_virus.fa eg2.sorted.bam | bcftools view -bvcg – > eg2.raw.bcf
5，看突變位點
bcftools view eg2.raw.bcf

七、討論

如何讓bowtie快起來：
如果bowtie在你的機器上執行起來很慢，那麼你可以試試以下的一些辦法來讓它跑得快一些：
1，儘可能的使用64位bowtie。很顯然，64位運算會比32位運算更快。所以最好使用支援64位運算的計算機來執行64位的bowtie。如果你是從原檔案開始編譯程式，在g++編譯時，你需要傳遞-m64引數，你也可以在make的時候加入這一資訊，比如說傳遞BITS=64給make，具體的：make BITS=64 bowtie。想知道你自己運算的是什麼版本的bowtie，你可以執行bowtie –version
2，如果你的計算機有多個CPU或者CPU核心，那麼請使用-p引數。-p引數會讓bowtie進入多執行緒模式。每一個執行緒都會使用單獨的CPU或者CPU核心。這種並行的運算模式也會大大加快運算速度。
3，如果你的報告檔案中每條短序列都有太多的匹配位點（超過10）那麼你可以試著重新使用bowtie-build加上 –offrate引數，如bowtie-build –offrate 4。-o/–offrate預設值為5，每下降1，比對速度會增加1倍，但是系統消耗（硬碟空間和記憶體）也會加倍。如果你的系統 配置太低，比如記憶體不足4GB，那麼建議你在bowtie的時候使用–offrate引數。與上一條相反的，你需要加大offrate的值。bowtie –offrate 6. 其預設值為5。每增加1，記憶體空間的要求下降，這樣會減少讀取硬碟當�%AL��擬記憶體的次數，速度反而會有所上升。
4，-n模式與-v模式。
預設的，bowtie採用了和Maq一樣的質量控制策略，設定-n 2 -l 28 -e 70。總的來說，比對模式分為兩種，一種是 -n 模式， 一種是 -v 模式，而且這兩種模式是不能同時使用的。bowtie預設使用-n模式。
-n模式引數： -n N -l L -e E
其中Ｎ，Ｌ，Ｅ都為整數。-n N 代表在高保真區內錯配不能超過Ｎ個，可以是0?3，一般的設定為2。-l L代表序列高保真區的長度，最短不能少於5，對於短序列長度為32的，設定為28就很不錯。-e E代表在錯配位點Phred quality值不能超過E，預設值為40。Phred quality值的計算式為：-10 log(P,base(10))
Phred Quality值 錯配可能性 正配可能性
10 1/10 90％
20 1/100 99％
30 1/1000 99.9％
40 1/10000 99.99%
50 1/100000 99.999%
而-v模式的引數相對較少。
-v模式引數：-v V
其中V為整數。-v V代表全長錯配不能超過Ｖ個，可以是0?3。這時，不考慮是否高保真區，也不考慮Phred quality值。
5，–best 與–strata
–best 引數代表報告檔案中，每個短序列的匹配結果將按匹配質量由高到低排序。–strata引數必須與–best引數一起使用，其作用是隻報告質量最高的那部 分。所謂質量高低，其實就是指錯配的鹼基數，如果指定了-l L引數，那就是在高保真區內的錯配數，否則就是全序列的錯配數。如果你還指定了 -M X的話，那就會在質量最高的當中，隨機選擇X個來報告。也就是說，當我們指定了-M 1 –best –strata的話，那就只報告1個最好的。
對於輸入，-q是指輸入的檔案為FASTQ(副檔名通常為.fq或者.fastq)格式；-f是指輸入檔案為FASTA(副檔名通常為.fa, .mfa或者.fna)格式；-c是指在命令列直接輸入要比對的序列。

參考資料：
官方使用說明手冊 http://computing.bio.cam.ac.uk/local/doc/bowtie2.html
bowtie:短序列比對的新工具http://blog.sina.com.cn/s/blog_5ecfd9d90100ukqq.html
Lemon的部落格：http://blog.163.com/[email protected]/blog/static/3530820120137132216950/

附：
[Bowtie 2] is often the first step in pipelines for comparative genomics,including for variation calling, ChIP-seq, RNA-seq, BS-seq. [Bowtie 2] and [Bowtie] (also called “[Bowtie 1]” here) are also tightly integrated into some
tools, including [TopHat]: a fast splice junction mapper for RNA-seq reads,[Cufflinks]: a tool for transcriptome assembly and isoform quantitiation from RNA-seq reads, [Crossbow]: a cloud-enabled software tool for analyzing

reseuqncing data, and [Myrna]: a cloud-enabled software tool for aligning RNA-seq reads and measuring differential gene expression.