【6.2.1】BBTools/BBMap Suite 的使用

更多参考:https://blog.csdn.net/msw521sg/article/details/52565985

BBMap/BBTools是一种用于DNA和RNA测序reads的拼接感知全局的比对工具。它可以处理的reads包括,Illumina,454,Sanger,Ion Torrent,Pac Bio和Nanopore。 BBMap快速且极其准确,特别是对于高度突变的基因组或长插入读数,甚至超过100kbp长的全基因缺失。它对基因组大小或contigs数量没有上限。并且已经成功地用于绘制到具有超过2亿个contigs的85gb的土壤宏基因组。此外,与其他比对工具相比,索引阶段非常快。 BBMap有很多选项,在shell脚本中描述。它可以输出许多不同的统计文件,例如经验读取质量直方图,插入大小分布和基因组覆盖,有或没有生成sam文件。因此,它可用于文库的质量控制和测序运行,或评估新的测序平台。衍生程序BBSplit也可用于分类或精炼宏基因读取。

那问题来了,现在的工具那么多为什么会推荐还有写这个bbmap、bbTools的教程?那当然是因为该工具有一些它自身的特点:

  1. BBTools是用纯Java编写的,可以在任何平台(PC,Mac,UNIX)上运行,除了安装Java之外没有依赖项(为Java 6及更高版本编译)。所有工具都是高效且多线程的。
  2. BBTools可以处理常见的排序文件格式,如fastq,fasta,sam,scarf,fasta + qual,压缩或raw,具有自动检测质量编码和交错。
  3. BBTool套件包括对高通量测序数据进行质量控制,修剪(trim),纠错,组装和比对的程序。
  4. BBTools是开源的,可以无限制免费使用,不用钱当然好!
  5. BBTools无需安装。从BBTools SF站点下载tar存档后,您需要做的就是解压,很动心是不是?

一、安装

方法一:

conda install -y bbmap

方法二:

##下载
wget https://excellmedia.dl.sourceforge.net/project/bbmap/BBMap_38.16.tar.gz
##解压
tar -xvfz BBMap_38.16.tar.gz

注意:BBMap(和所有相关工具)shellcripts将尝试自动检测内存,但可能会失败(导致jvm无法启动或内存不足),具体取决于系统配置。这可以通过将-XmxNNg标志添加到shellscript的参数列表(根据计算机的物理内存进行调整,不要使用超过总RAM的80%)来覆盖,并将其传递给java。

下载地址: https://sourceforge.net/projects/bbmap/

二、使用

BBTools程序接受经过gzip压缩的fastq文件,并能够以多种数据格式写入输出。同时支持以fq 或者 fastq结尾的测序文件。

2.1 统计序列

stats (stats.sh)

用于统计reads或者assembly基本的信息,例如 scaffold count, N50, L50, GC content, gap percent。

#用一个fastq文件作为例子,这个工具也可以用于fasta文件
stats.sh in=A1.1.fastq

#结果
Main genome scaffold total:             250000
Main genome contig total:               250000
Main genome scaffold sequence total:    22.500 MB
Main genome contig sequence total:      22.500 MB       0.000% gap
Main genome scaffold N/L50:             250000/90
Main genome contig N/L50:               249990/90
Main genome scaffold N/L90:             250000/90
Main genome contig N/L90:               249990/90
Max scaffold length:                    90
Max contig length:                      90
Number of scaffolds > 50 KB:            0
% main genome in scaffolds > 50 KB:     0.00%


Minimum         Number          Number          Total           Total           Scaffold
Scaffold        of              of              Scaffold        Contig          Contig
Length          Scaffolds       Contigs         Length          Length          Coverage
--------        --------------  --------------  --------------  --------------  --------
    All                250,000         250,000      22,500,000      22,499,990   100.00%
     50                250,000         250,000      22,500,000      22,499,990   100.00%

如果你有多个输入文件,可以使用statswrapper.sh工具,对多个数据同时进行运算:

statswrapper.sh *.fastq

n_scaffolds     n_contigs       scaf_bp contig_bp       gap_pct scaf_N50        scaf_L50        ctg_N50 ctg_L50 scaf_N90        scaf_L90        ctg_N90 ctg_L90 scaf_max        ctg_max scaf_n_gt50K    scaf_pct_gt50K  gc_avg  gc_std filename
250000  250000  22500000        22499990        0.000   250000  90      249990  90      250000  90      249990  90      90      90      0       0.000   0.48428 0.05525 /scratch/pawsey0149/hhu/bio_trainee/tool_bbmap/data/A1.1.fastq
250000  250000  22500000        22499464        0.002   250000  90      249703  90      250000  90      249703  90      90      90      0       0.000   0.48673 0.05146 /scratch/pawsey0149/hhu/bio_trainee/tool_bbmap/data/A1.2.fastq

该工具在统计多个fastq文件的基础信息时非常方便

pileup (pileup.sh)

用于计算scaffold的覆盖率,需要没有sort过的bam 或者 sam 的格式的文件作为input。

 pileup.sh in=test.sam out=test.out

 #屏幕输出 一些基本的比对信息,还有比对使用的ref的信息

Reads:                                  66234
Mapped reads:                           57247
Mapped bases:                           8587050
Ref scaffolds:                          9
Ref bases:                              41030279

Percent mapped:                         86.431
Percent proper pairs:                   0.000
Average coverage:                       0.209
Standard deviation:                     0.583
Percent scaffolds with any coverage:    100.00
Percent of reference bases covered:     15.45

#test.out里面的输出,每一个染色体上比对的覆盖率等信息。
#ID     Avg_fold        Length  Ref_GC  Covered_percent Covered_bases   Plus_reads      Minus_reads     Read_GC Median_fold     Std_Dev
1       0.2132  7978604 0.0000  15.6789 1250959 5765    5577    0.5149  0       0.59
2       0.2127  8319966 0.0000  15.6720 1303906 5940    5860    0.5195  0       0.59
3       0.2031  6606598 0.0000  15.0522 994440  4417    4531    0.5205  0       0.58
4       0.2047  5546968 0.0000  14.9938 831704  3795    3775    0.5218  0       0.58
5       0.2126  4490059 0.0000  15.6697 703581  3183    3182    0.5188  0       0.59
6       0.2124  4133993 0.0000  15.8143 653760  2973    2882    0.5202  0       0.58
7       0.2068  3415785 0.0000  15.3674 524917  2315    2394    0.5139  0       0.58
supercont8.8    0.1830  535760  0.0000  14.1227 75664   325     329     0.5165  0       0.52
mit     0.2357  2546    0.0000  23.5664 600     2       2       0.6050  0       0.42

readlength (readlength.sh)

统计read的平均长度,总长度等数据信息

readlength.sh in=A1.1.fastq in2=A1.2.fastq

#Reads: 500000
#Bases: 45000000
#Max:   90
#Min:   90
#Avg:   90.0
#Median:        90
#Mode:  90
#Std_Dev:       0.0
#Read Length Histogram:
#Length reads   pct_reads       cum_reads       cum_pct_reads   bases   pct_bases       cum_bases       cum_pct_bases
90      500000  100.000%        500000  100.000%        45000000        100.000%        45000000        100.000%

kmercountexact (kmercountexact.sh)

计算文件中unique kmers的数量。生成kmer频率直方图和基因组大小估计。

kmercountexact.sh in=A1.1.fastq in2=A1.2.fastq out=kmer_test.out khist=kmer.khist peaks=kmer.peak

##屏幕输出
Input:                          500000 reads            45000000 bases.

For K=31
Unique Kmers:                   3326507
Average Kmer Count:             9.016
Estimated Kmer Depth:           9.015
Estimated Read Depth:           13.526

##新生成三个文件
kmer_test.out ##运行输出,所以kmer的fa集合
kmer.khist ##成kmer频率直方图列表
kmer.peak ##kmer统计size等峰值的统计信息

bbmask (bbmask.sh)

Masks低复杂度或包含重复kmers的序列,或由比对到的reads所覆盖的序列。默认情况下的参数:a window=80 and entropy=0.75

输入可以是stdin或fasta或fastq文件,sam也是可选的文件格式。

### 其中一个例子
bbmask.sh in=A1.1.fastq out=test.mark

### 屏幕输出

Masking low-entropy (to disable, set 'mle=f')
Low Complexity Masking Time:    0.479 seconds.
Ref Bases:                  22500000    46.97m bases/sec
Low Complexity Bases:             88

Converting masked bases to N
Done Masking
###这里就是将fastq文件进行简单重复的屏蔽,将其标记为N,这里可以看到有88个被标记为N

2.2 序列处理

BBMap Read Merger

合并双末端(PE)reads,在预期中这些reads有重叠的位置(例如16S扩增子测序)。

通过重叠检测将配对的reads合并为单个read。如果具有足够的覆盖率,还可以通过kmer扩展合并非重叠的reads。 Kmer模式需要更多内存,应与bbmerge-auto.sh脚本一起使用。

Usage for interleaved files (交错的文件):    bbmerge.sh in=<reads> out=<merged reads> outu=<unmerged reads>
Usage for paired files (有双末端的文件):         bbmerge.sh in1=<read1> in2=<read2> out=<merged reads> outu1=<unmerged1> outu2=<unmerged2>

BBMap bbsplit

如何使用多个reference对我的样本中的reads进行多次比对?

这是一个很常见的问题,bbsplit是一个很好的工具去处理该问题。它非常灵活,可以处理多个数据序列。

同时比对多个基因组。此功能可用于将reads(可用于去污染)分类到“不同基因组”特定池中。

基本用法

To index:     bbsplit.sh build=<1> ref_x=<reference fasta> ref_y=<another reference fasta>
To map:       bbsplit.sh build=<1> in=<reads> out_x=<output file> out_y=<another output file>

然后这两个命令可以并在一起:

bbsplit.sh ref=x.fa,y.fa in=reads.fq basename=o%.fq

# this is equivalent to

bbsplit.sh build=1 in=reads.fq ref_x=x.fa ref_y=y.fa out_x=ox.fq out_y=oy.fq

默认情况下,配对读取将产生交错输出,但您可以使用#符号生成双输出文件。例如,basename = o%_#。fq将生成ox_1.fq,ox_2.fq,oy_1.fq和oy_2.fq。

我们现在将通过一个简单的例子来使用下面的bbsplit。

#工具下载
conda install entrez-direct

# 数据下载,获得来自NCBI的三种细菌基因组的fasta格式化序列。

efetch -db nucleotide -id NC_003454 -format fasta > NC_003454.fa
efetch -db nucleotide -id NC_010468 -format fasta > NC_010468.fa
efetch -db nucleotide -id NC_013520 -format fasta > NC_013520.fa

# 下载细菌的reads

wget http://data.biostarhandbook.com/reads/mixed-bacterial-reads.fq.gz

#
尽管我们在这用到的是细菌的reads,但这个方法通常适用于任何混合数据。我们可以使用一个或多个基因组的任意组合来将reads分类进每个对应的基因组中。

bbsplit.sh in=mixed-bacterial-reads.fq.gz ref=NC_010468.fa,NC_013520.fa,NC_003454.fa basename=out_%_#.fq.gz outu=rest.fq.gz minratio=0.76 minhits=1

bbsplit运行后应生成以下文件。

ls -sh1tr

2.4M rest.fq.gz
608K out_NC_013520_1.fq.gz
608K out_NC_003454_1.fq.gz
608K out_NC_010468_1.fq.gz

mixed-bacterial-reads.fq.gz文件现已被分成三个文件,其中包含我们感兴趣的reads。剩余的reads放在rest.fq.gz中

BBMap repair

配对末端(PE)序列是一种非常有用的工具,可提供有关被测序片段的空间信息。 PE的默认输出将读取对(指定为R1和R2)放在两个单独的文件中。当您扫描这些文件是否存在adapter序列和其他污染物时(相当于trimming,做质控的时候),建议将PE数据文件一起处理。 PE识别修剪/扫描程序使两个reads文件的reads顺序保一致。文件内部的读取顺序很重要,因为大多数NGS的mapping 工具都希望两个文件中的reads顺序相同(并且可能无法检查)。

如果文件被错误地trim或者某些文件被破坏,PE中的读取顺序可能会受到干扰。在最好的情况下,当ni使用此类文件作为NGS mapping 工具的输入时,这将产生错误。 (注意:在最坏的情况下,mapping可能会按原样处理数据。这可能会导致大量不一致的比对产生)。

使用BBMap工具中的repair.sh可以“重新同步保持一致”PE文件。这是一个举例说明的例子:

首先,人为先创造出一个不配对的fastq双末端的文件。

#下次会继续详细介绍这个工具,这里先使用
reformat.sh in=A1.2.fastq samplerate=0.9 out=A1.2.bad.fastq

这样子,就会产生一组不对称的双末端fastq file,一个是A1.1.fastq 另一个是A1.2.bad.fastq

接下来我们准备用repair.sh来修复’损坏的’reads配对并重新同步PE reads 以生成“固定”reads文件。 repair.sh的输出是交错的PE的reads(read1_R1,read1_R2,read2_R1,read2_R2等)。 “破碎”的reads将被存储于A1.bad.fastq中:

repair.sh in1=A1.1.fastq in2=A1.2.bad.fastq out=A1_fixed.fastq outsingle=A1.bad.fastq

我们现在可以通过reformat.sh 这个工具,将固定reads文件A1_fixed.fastq拆分为单独的R1 / R2 reads文件:

reformat.sh in=A1_fixed.fastq out1=A1.1.fixed.fastq out2=SRR1972739_2_fixed.fastq addcolon

简单来说,这个小工具是处理两个不一致的PE fastq文件的小神器。

2.3 BBMap suite - Simulation of data (数据模拟)

除了之前介绍的工具,BBMap还可以用来, 从参考基因组产生随机合reads。reads的名称表示其基因组来源。并可以允许准确定的insert size和合成突变的类型,大小和速率等具体内容。 Illumina和PacBio类型的reads都可以使用randomreads这个小工具来生成。

您可以指定PE 的reads,insert size的分布,reads长度(或长度范围)等。 Randomreads专门设计用于对突变进行出色的控制。你可以准确地或概率地指定每个reads所包含的snps,insertion,deletions和Ns的数量;这些突变的长度可以进行单独的定制,质量值可以交替设置,以允许根据质量生成错误。所有reads都使用其基因组来源进行注释。

请记住,50%的reads将来自正链,50%来自负链。因此,reads将完全匹配参考基因组,或者它的反向补链将完美匹配。

randomreads.sh ref=<file> out=<file> length=<number> reads=<number>

介绍那么多,可能大家还是会有点疑惑究竟,这个工具有啥用?废话不多说,通过几行代码带你快速了解该工具的用处:

模拟高通量测序Illumina数据集。

1.首先先下载序列号为NC_010468 fasta格式的文件

efetch -db=nuccore -format=fasta -id=NC_010468 > NC_010468.fa

我们将通过以下命令简化NC_010468.fa中的fasta标头,以便新标头读取> NC_010468_E_coli

sed '1 s/^.*$/\>NC_010468_E_coli/' NC_010468.fa > NC_010468_new.fa

我们将使用此参考序列文件生成的合成illumina配对末端数据集(2 x 300 bp)

randomreads.sh ref=NC_010468_new.fa out=NC_010468_synth.fq.gz length=300 paired=t mininsert=100 maxinsert=400 reads=1000000

#默认情况下,生成的reads采用交错格式。它们可以通过`reformat.sh`进一步分成单独的read文件

reformat.sh in=NC_010468_synth.fq.gz out1=NC_010468_synth_R1.fq.gz out2=NC_010468_synth_R2.fq.gz addcolon

2.使用randomreads生成PacBio格式数据集(使用PacBio错误模型)。 模拟一个包含10,000个读数的PacBio数据集,其最小读取长度为1kb,最大读取长度为5kb。

randomreads.sh ref=NC_010468_new.fa out=NC_010468_synth_pacbio.fq.gz minlength=1000 maxlength=5000 reads=10000 pacbio=t

3.生成更复杂/有趣的合成数据集,可以模拟宏基因组。 下载数据

efetch -db=nuccore -format=fasta -id=NC_009614 > NC_009614.fa
efetch -db=nuccore -format=fasta -id=NC_013520 > NC_013520.fa

metagenome将为每个scaffold分配一个随机指数变量,该变量决定从该scaffold生成read的概率。因此,如果将20个细菌基因组连接在一起,则可以运行随机read并获得类似宏基因组的分布。当使用转录组参考时,它也可以用于RNA-seq。

覆盖率取决于每个参考序列水平,因此如果细菌组装具有多个contigs,您可能需要先将它们与fuse.sh粘在一起,然后再与其他参考物连接。

首先将ref 合成单一的文件:

cat NC_010468.fa NC_013520.fa NC_009614.fa AF086833.fa > metag.fa

然后我们将使用它作为宏基因组生成的输入,然后使用reformat.sh将读取分成两个文件。

randomreads.sh ref=metag.fa out=stdout length=300 paired=t mininsert=100 maxinsert=400 reads=5000000 metagenome=t coverage=3 | reformat.sh in=stdin out1=metag_R1.fq.gz out2=metag_R2.fq.gz interleaved addslash

2.4 重新格式化和过滤序列数据

Reformat.sh这个工具包前面一直有用到,接下来会详细讲解一下他的使用。

功能: 格式转换/操作/分析序列数据。

重新格式化reads以更改ASCII质量编码,交错,文件格式或压缩格式。可选择执行其他功能,如质量修整,子集和子采样。支持fastq,fasta,fasta,bam,gzip 等得格式。

reformat.sh in=<file> in2=<file2> out=<outfile> out2=<outfile2>
#in2 and out2 这两个options是PE reads的文件才会用到的(双末端)

具体应用

1.切除一部分的reads文件

使用reformat.sh可以轻松地从文件中产生随机获取的reads。

#随机抽取300个reads
reformat.sh in=SRR1972739_1.fastq out=SRR1972739_1_sample.fastq sample=300


#您可以通过执行以总reads的百分比进行采样(以下示例中为10%):
reformat.sh in=SRR1972739_1.fastq out=SRR1972739_1_10perc.fastq samplerate=0.1

2.验证配对端数据的reads是正确配对的。

reformat.sh in1=r1.fq.gz in2=r2.fq.gz vpair

3.按reads的长度过滤SAM / BAM文件

reformat.sh in=x.sam out=y.sam minlength=50 maxlength=200

4.过滤/检测q SAM / BAM文件中的拼接reads 您可以将maxdellen(在下面的示例中为50)设置为你认为最小的任何长度的eletion事件,以表示拼接,这取决于你的物种的类型。

reformat.sh in=mapped.bam out=filtered.bam maxdellen=50

5.从SAM / BAM文件中获取比对不上的reads

reformat.sh in=all.sam out=unmapped.fq.gz unmappedonly
repair.sh in=unmapped.fq.gz  out=r1.fq.gz out2=r2.fq.gz outs=singleton.fq

2.5 删除重复的序列数据集

从一组文件中删除重复的读取/序列。 类似picard tool的Markduplicates的功能。

dedupe.sh in=<file or stdin> out=<file or stdout>

在运行完成后,还会生成一组全面的统计信息。它们通常被写入\但可以很容易地捕获到日志文件中(例如下面的例子)

Input:                          15000 reads             1515000 bases.
Duplicates:                     7445 reads (49.63%)     751945 bases (49.63%)           0 collisions.
Containments:                   0 reads (0.00%)         0 bases (0.00%)         113772 collisions.
Result:                         7555 reads (50.37%)     763055 bases (50.37%)

四、我的案例

参考资料

个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn

Sam avatar
About Sam
专注生物信息 专注转化医学