【4.1.1】bowtie2

一、简介

Bowtie是一个超级快速的,较为节省内存的短序列拼接至模板基因组的工具。它在拼接35碱基长度的序列时,可以达到每小时2.5亿次的拼接速度。 Bowtie并不是一个简单的拼接工具,它不同于Blast等。它适合的工作是将小序列比对至大基因组上去。它最长能读取1024个碱基的片段。模板最小尺寸不能小于1024碱基,而短序列最长而不能超过1024碱基。换言之,bowtie非常适合下一代测序技术。 在 使用bowtie前,需要使用bowtie-build来构建比对模板。

1.1 转录组还是基因组?

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也可以调用它

1.2 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. Ns) in the reference. Bowtie 1 does not.
  6. Bowtie 2不能比对colorspace reads.

1.3 选不选bowtie

  • Bowtie 2 在比对大的基因组的时候,工作效果更好一些
  • 如果你是比对两个特别大的基因组,可以考虑MUMmer
  • 如果你比对的是两个特别相近且比较小的基因组,可以用NUCmer, BLAT, or BLAST
  • 但是Bowtie 2不能比对colorspace reads.

1.4 Bowtie设计思路

  1. 短序列在基因组上至少有一处最适匹配,
  2. 大部分的短序列的质量是比较高
  3. 短序列在基因组上最适匹配的位置最好只有一处。这些标准基本上和RNA-seq, ChIP-seq以及其它一些正在兴起的测序技术或者再测序技术的要求一致。

1.5 使用tophat的基本步骤是:

使用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下载地址

cd /data/software/

wget -c https://sourceforge.net/projects/bowtie-bio/files/bowtie/1.0.1/bowtie-1.0.1-linux-x86_64.zip

unzip bowtie-1.0.1-linux-x86_64.zip

cd /data/software/tools/bowtie-1.0.1

vim /etc/profile
export PATH=/data/software/tools/bowtie-1.0.1:$PATH

bowtie2下载地址

https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.5.1/

cd /data/software/tools/
wget -c https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.5.1/bowtie2-2.3.5.1-linux-x86_64.zip

/data/software/tools/bowtie2-2.3.5.1-linux-x86_64

unzip bowtie2-2.3.5.1-linux-x86_64.zip

vim /etc/profile
export PATH=/data/software/tools/bowtie2-2.3.5.1-linux-x86_64:$PATH

三、使用案例

3.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文件的时候,可以用逗号分开处理

3.2 针对unpaired reads的比对

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

3.3 针对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就是你的选择。

3.4 局部比对的例子

用[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

3.5 其他参数的说明:

-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 怎么设置种子,详见官网说明书

…… 其他的参数具体看官网说明书

四、讨论

4.1 如何让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,L,E都为整数。-n 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代表全长错配不能超过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是指在命令行直接输入要比对的序列。

4.2 使用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

五、报错

5.1 bowtie2-align died with signal 6

sam@sam-Precision-WorkStation-T7500[index] /sam/bowtie2/bowtie2 -x scaffold -1 /sam/idba_ud/myproject/reads2/total_clean.1.fq -2 /sam/idba_ud/myproject/reads2/total_clean.2.fq -S eg2.sam
Error, fewer reads in file specified with -2 than in file specified with -1
terminate called after throwing an instance of 'int'
bowtie2-align died with signal 6 (ABRT)
I use trimmomatic to trim my PE reads data ,so get 4 output,I cat the output together:
cat output_forward_paired.fq output_forward_unpaired.fq >total_clean.1.fq
cat output_reverse_paired.fq output_reverse_unpaired.fq >total_clean.2.fq

这样的错误源于我对官网错误的理解,官网原文如下:

官网:http://computing.bio.cam.ac.uk/local/doc/bowtie2.html Mixed mode: paired where possible, unpaired otherwise If Bowtie 2 cannot find a paired-end alignment for a pair, by default it will go on to look for unpaired alignments for the constituent mates. This is called “mixed mode.” To disable mixed mode, set the –no-mixed option. Bowtie 2 runs a little faster in –no-mixed mode, but will only consider alignment status of pairs per se, not individual mates. I.e. in SAM records for pairs, the optional YM:i flag will not appear and MAPQ values calculated for mates may be underestimates in some cases.

我误以为这个bowtie既可以处理pair的,也可以同时处理unpair的,结果出现上面的问题。 实际上应该是bowtie2处理的 paired reads 在的两个文件,必须有同样多的 reads数,并且都是一一对应的;如果一对不能都比上 会找其中一条来比对上;一个如果不能比对到基因组,那么使用 –no-mixed 参数,则不会对与之匹配的reads进行比对了。总之核心思想就是必须处理一对,如果一对中有一条未能匹配上,你可以选择是否保留这一对的结果,–no-mixed就是你的选择。

那如何同时处理pair和unpair?

先用paired 去比对

/sam/bowtie2/bowtie2 -p 12 -x scaffold -1 /sam/idba_ud/myproject/reads2/output_forward_paired.fq -2 /sam/idba_ud/myproject/reads2/output_reverse_paired.fq -S eg2.sam

再用unpair去比对

/sam/bowtie2/bowtie2 -p 12 -x scaffold -U /sam/idba_ud/myproject/reads2/unpair.fq -S eg1.sam

接下来就注意了,如果直接去cat sam文件的话,因为sam文件有开头的标志,所以后面的生成bam文件就会出错,这里就有三种思路可以来解决:

cat sam文件(但得去掉第二个文件sam文件标志性的开头)
samtools merge  针对排序后的文件
samtools cat 针对bam文件

这里面我犯的错误如下:

cat eg1.sam eg2.sam >eg.sam
samtools cat eg1.sam eg2.sam >eg.sam
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
zsh: segmentation fault  samtools cat eg1.sam eg2.sam > eg.sam

samtools cat eg1.sam eg2.sam >eg.sam
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
zsh: segmentation fault  samtools cat eg1.sam eg2.sam > eg.sam

我采用的第三种方式

samtools view -bS eg1.sam > eg1.bam
samtools view -bS eg2.sam > eg2.bam  
samtools cat eg1.bam eg2.bam >eg.bam

ps:

  • 问题能真相大白,首先感谢QQ好友陈连福和哈皮耐心的讲解
  • 同时bowtie2创始人的邮箱:langmea@cs.jhu.edu和sourceforge上也有论坛
  • sourceforge论坛回复地址:https://sourceforge.net/p/bowtie-bio/bugs/305/

参考资料

  • 官方使用说明手册 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 reseuqncing data, and [Myrna]: a cloud-enabled software tool for aligning RNA-seq reads and measuring differential gene expression.

这里是一个广告位,,感兴趣的都可以发邮件聊聊:tiehan@sina.cn
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn