【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的教程？那当然是因为该工具有一些它自身的特点：
- BBTools是用纯Java编写的，可以在任何平台（PC，Mac，UNIX）上运行，除了安装Java之外没有依赖项（为Java 6及更高版本编译）。所有工具都是高效且多线程的。
- BBTools可以处理常见的排序文件格式，如fastq，fasta，sam，scarf，fasta + qual，压缩或raw，具有自动检测质量编码和交错。
- BBTool套件包括对高通量测序数据进行质量控制，修剪(trim)，纠错，组装和比对的程序。
- BBTools是开源的，可以无限制免费使用，不用钱当然好！
- 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%)
四、我的案例
参考资料
