简单介绍不会编程的小白通常会烦恼怎样去处理一些简单的数据统计,不要紧其实日常我们用到百分之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的教程?那当然是因为该工具有一些它自身的特点:
BBTools是用纯Java编写的,可以在任何平台(PC,Mac,UNIX)上运行,除了安装Java之外没有依赖项(为Java 6及更高版本编译)。所有工具都是高效且多线程的。
BBTools可以处理常见的排序文件格式,如fastq,fasta,sam,scarf,fasta + qual,压缩或raw,具有自动检测质量编码和交错。
BBTool套件包括对高通量测序数据进行质量控制,修剪(trim),纠错,组装和比对的程序。
BBTools是开源的,可以无限制免费使用,不用钱当然好!
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: 66234Mapped reads: 57247Mapped bases: 8587050Ref scaffolds: 9Ref bases: 41030279Percent mapped: 86.431Percent proper pairs: 0.000Average coverage: 0.209Standard deviation: 0.583Percent scaffolds with any coverage: 100.00Percent 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_Dev1 0.2132 7978604 0.0000 15.6789 1250959 5765 5577 0.5149 0 0.592 0.2127 8319966 0.0000 15.6720 1303906 5940 5860 0.5195 0 0.593 0.2031 6606598 0.0000 15.0522 994440 4417 4531 0.5205 0 0.584 0.2047 5546968 0.0000 14.9938 831704 3795 3775 0.5218 0 0.585 0.2126 4490059 0.0000 15.6697 703581 3183 3182 0.5188 0 0.596 0.2124 4133993 0.0000 15.8143 653760 2973 2882 0.5202 0 0.587 0.2068 3415785 0.0000 15.3674 524917 2315 2394 0.5139 0 0.58supercont8.8 0.1830 535760 0.0000 14.1227 75664 325 329 0.5165 0 0.52mit 0.2357 2546 0.0000 23.5664 600 2 2 0.6050 0 0.42readlength (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_bases90 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=31Unique Kmers: 3326507Average Kmer Count: 9.016Estimated Kmer Depth: 9.015Estimated 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
作者:lakeseafly
链接:https://www.jianshu.com/p/175c3282a61c



