转载：http://sega97.blog.163.com/blog/static/7604122014229105912100/
http://www.plob.org/2012/09/02/4540.html
转录组还是基因组？

map常用的工具有bowtie/bowtie2, BWA,SOAP1/SOAP2等。这个问题又会被分成两个问题，是基因组测序（DNA-seq）还是转录组测序(mRNA-seq)。其中的区别是对于真核生物而言，mRNA序列与DNA序列并不完全相同，在经历了后剪切之后，成熟的mRNA可能是原基因的一部分，甚至顺序及个别碱基会产生变化。如果是mRNA测序，那map工作就会在DNA测序map的基础上再多一步，map到转录组上去。所以最为流行的做法是，使用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，

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

 

其他参数的说明：

    -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   同理同上，在我之前的博客中有记录

    --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/zhoulili1987619@126/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