生信小工具专题:BBTools/BBMap Suite 的使用 (1)

不会编程的小白通常会烦恼怎样去处理一些简单的数据统计,不要紧其实日常我们用到百分之80以上的数据统计分析,都已经有人帮我们将轮子造好了,只要我们学会怎样利用好这些工具,很多事情都可以事半功倍。这周开始会给大家介绍一些常用的生信工具,并且提供一些简单的测试数据供大家学习玩耍。第一篇文章首先讲BBMap。

简单介绍

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存档后,您需要做的就是解压,很动心是不是?

安装与注意事项

经过简单的介绍后,大家相信对这一款工具也有一定的了解,下面开始安装该工具。这里的安装介绍都是基于Linux系统,Windows系统的运行方式有点不同,在这里没有介绍。

conda 大法好,简单一个command帮你搞定安装,如果你还不知道什么是conda,可以先阅读一下大神的博客 (http://www.bio-info-trainee.com/1906.html),如果你还是没有弄懂,或者你想了解更多关于文件安装之类的话题,那么我强力推荐你去观看由我们大号生信技能树所推出的一个课程:

conda install -y bbmap

当然如果你不会使用conda,也没问题,直接去到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。

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

牛试小刀

在按照工具类别对BBmap/BBtools进行逐步介绍时,先介绍几款其在统计数据中简单但有好用的小工具。

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

今天就简单介绍到这,有条件的同学可以用我提供的测试数据进行理解与学习。下一个推文会按照BBmap工具功能进行分类并详细介绍相应的工具功能。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 157,012评论 4 359
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 66,589评论 1 290
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 106,819评论 0 237
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,652评论 0 202
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 51,954评论 3 285
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,381评论 1 210
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,687评论 2 310
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,404评论 0 194
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,082评论 1 238
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,355评论 2 241
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,880评论 1 255
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,249评论 2 250
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,864评论 3 232
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,007评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,760评论 0 192
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,394评论 2 269
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,281评论 2 259

推荐阅读更多精彩内容