BBtoolsはアメリカのJGIが提供している多機能なNGS向けの解析ツール。2014年にオープンソース化されたらしい。論文は現在準備中とある。アライメントのBBmapや、オーバーラップがないペアリードをマージするBBMerge、エラーコレクションしたfastqを出力するBBnormなど特徴的なコマンドを備える。ここでは以下のコマンドを紹介する。
- BBMap - リファンレンスへのアライメント。
- BBMerge - ペアリードのマージ。k-mer配列を付け足し延長することも可能。
- BBnorm - カバレッジのノーマライズ。
- Dedupe - duplicateしたcontigの除去。
其の2。
公式ページ
http://jgi.doe.gov/data-and-tools/bbtools/
ユーザーガイド
http://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/
ダウンロード
https://sourceforge.net/projects/bbmap/
またはbrewでも導入できる。
brew install BBtools
BBtoolsのコマンドは指定がない限りコア数をあるだけ使う。
ラン
BBMap (リンク)-リードのリファンレンスへのアライメント
BBMapは、k-merに分解してk-merサイズのseedでアライメントと伸長を行うアライナー(デファルトではk=13)。 ショートリードのアライメントだけでなく、ロングリードのアライメントも可能とされ、Illumina,、454、Sanger、Ion Torrent、Pac Bio、Nanoporeのリードのアライメント対応がマニュアルに記載されている。
・index作成とペアリードのアライメント(index=>マッピング=>アライメント)
bbmap.sh in1=reads1.fq in2=reads2.fq out=mapped.sam ref=reference.fa
indexディレクトリが作られ、リードがリファレンスにアライメントされる。defaultだとマシンの全コアが計算に使用される(明示するなら"t=")。
bbmap.shは700bp以下のリードしか扱えない。pacbioのアライメントは6-kbまでのリードに対応した専用コマンドを使う。
・ロングリードのアライメント。
mapPacBio.sh in=reads.fq out=mapped.sam ref=reference.fa
6-kb以上のリードは自動で断片化されリネームされてアライメントされる。ナノポアのデータもこのコマンドを使うことが推奨されている。
・様々な統計データを出力する。
bbmap.sh in1=reads1.fq in2=reads2.fq bhist=bhist.txt qhist=qhist.txt aqhist=aqhist.txt lhist=lhist.txt ihist=ihist.txt ehist=ehist.txt qahist=qahist.txt indelhist=indelhist.txt mhist=mhist.txt gchist=gchist.txt idhist=idhist.txt scafstats=scafstats.txt
標準出力には以下のようなデータがプリントされる。
------------------ Results ------------------
Genome: 1
Key Length: 13
Max Indel: 16000
Minimum Score Ratio: 0.56
Mapping Mode: normal
Reads Used: 2647748 (397162200 bases)
Mapping: 70.725 seconds.
Reads/sec: 37437.34
kBases/sec: 5615.60
Pairing data: pct reads num reads pct bases num bases
mated pairs: 99.9427% 1323116 99.9427% 396934800
bad pairs: 0.0178% 236 0.0178% 70800
insert size avg: 601.02
insert 25th %: 566.00
insert median: 600.00
insert 75th %: 634.00
insert std dev: 52.83
insert mode: 602
Read 1 data: pct reads num reads pct bases num bases
mapped: 99.9677% 1323447 99.9677% 198517050
unambiguous: 98.0062% 1297479 98.0062% 194621850
ambiguous: 1.9615% 25968 1.9615% 3895200
low-Q discards: 0.0000% 0 0.0000% 0
perfect best site: 38.6619% 511835 38.6619% 76775250
semiperfect site: 38.6630% 511849 38.6630% 76777350
rescued: 1.2452% 16485
Match Rate: NA NA 99.3118% 197215279
Error Rate: 61.3256% 811612 0.6882% 1366602
Sub Rate: 60.9806% 807046 0.6527% 1296069
Del Rate: 0.5644% 7469 0.0326% 64831
Ins Rate: 0.3730% 4937 0.0029% 5702
N Rate: 0.0000% 0 0.0000% 0
Read 2 data: pct reads num reads pct bases num bases
mapped: 99.9702% 1323479 99.9702% 198521850
unambiguous: 97.9981% 1297371 97.9981% 194605650
ambiguous: 1.9721% 26108 1.9721% 3916200
low-Q discards: 0.0000% 0 0.0000% 0
perfect best site: 38.6839% 512126 38.6839% 76818900
semiperfect site: 38.6856% 512148 38.6856% 76822200
rescued: 1.1880% 15728
Match Rate: NA NA 99.3145% 197221255
Error Rate: 61.3046% 811353 0.6855% 1361264
Sub Rate: 60.9506% 806669 0.6523% 1295263
Del Rate: 0.5679% 7516 0.0306% 60669
Ins Rate: 0.3843% 5086 0.0027% 5332
N Rate: 0.0000% 0 0.0000% 0
Total time: 74.777 seconds.
さらに指定した分析データ(ヒストグラム)が出力される。
参考1
index作成時のメモリ使用量はリファンレス1塩基につき6 byteとされている。ヒトゲノムでは24GB必要(usemodulo flagをつけて12GB)。
・defaultのk=13の時のメモリ使用量は以下のコマンドで確認できる。
stats.sh in=reference.fa k=13
user$ stats.sh in=../testgenome.fasta k=13
A C G T N IUPAC Other GC GC_stdev
0.2629 0.2364 0.2372 0.2635 0.0000 0.0000 0.0000 0.4736 0.0241
Main genome scaffold total: 8
Main genome contig total: 8
Main genome scaffold sequence total: 3.957 MB
Main genome contig sequence total: 3.957 MB 0.000% gap
Main genome scaffold N/L50: 1/3.573 MB
Main genome contig N/L50: 1/3.573 MB
Main genome scaffold N/L90: 1/3.573 MB
Main genome contig N/L90: 1/3.573 MB
Max scaffold length: 3.573 MB
Max contig length: 3.573 MB
Number of scaffolds > 50 KB: 4
% main genome in scaffolds > 50 KB: 98.63%
Minimum Number Number Total Total Scaffold
Scaffold of of Scaffold Contig Contig
Length Scaffolds Contigs Length Length Coverage
-------- -------------- -------------- -------------- -------------- --------
All 8 8 3,956,956 3,956,956 100.00%
1 KB 8 8 3,956,956 3,956,956 100.00%
2.5 KB 6 6 3,952,233 3,952,233 100.00%
5 KB 6 6 3,952,233 3,952,233 100.00%
10 KB 5 5 3,947,019 3,947,019 100.00%
25 KB 5 5 3,947,019 3,947,019 100.00%
50 KB 4 4 3,902,676 3,902,676 100.00%
100 KB 4 4 3,902,676 3,902,676 100.00%
250 KB 1 1 3,573,470 3,573,470 100.00%
500 KB 1 1 3,573,470 3,573,470 100.00%
1 MB 1 1 3,573,470 3,573,470 100.00%
2.5 MB 1 1 3,573,470 3,573,470 100.00%
BBMap minimum memory estimate at k=13: -Xmx890m (at least 1000 MB physical RAM)
使用したバクテリアゲノムだと1GBの物理メモリが必要と算出された。
参考2
indel部位のスプリットアライメントに対応している。条件はmaxindelフラッグで規定されており、デフォルトではmaxindel=16000となっている。RNA seqでlarge intronが想定される場合などはmaxindel=100000などの大きな数値に指定することが推奨されている。
BBMerge (リンク)-ペアリードのマージ。
BBMegeはオーバーラップするペエアリードをマージして1つの長いリードを作るツール。
・ペアリードのマージ。
bbmerge.sh in1=reads1.fq in2=reads2.fq out=merged.fq outu1=unmerged1.fq outu2=unmerged2.fq ihist=ihist.txt
- in=null Primary input. 'in2' will specify a second file.
- out=<file> File for merged reads. 'out2' will specify a second file.
- outu=<file> File for unmerged reads. 'outu2' will specify a second file.
- ihist=<file> (hist) Insert length histogram output file.
インサートサイズのヒストグラムがihist.txtに出力される。マージを厳密に行うならpfilter=1をつける(“pfilter=1” is the strictest possible setting, which bans merges in which there are any mismatches in the overlap region. )。または、verystrict(Greatly decrease FP and merging rate.)をつける。オーサーはアセンブルに使うデータならverystrictをつけ、とにかくマージするならloose(Increase false positive rate and merging rate)をつけることを提案している。
・オーバーラップするペアリードのエラー修復。
bbmerge.sh in=interlace.fq out=corrected.fq ecco mix
- mix=f Output both the merged (or mergable) and unmerged reads in the same file (out=). Useful for ecco mode.
- ecco=f Error-correct the overlapping part, but don't merge.
ここでは入力するペアリードはインターレースとなっている。このコマンドを実行することでオーバーラップする領域で塩基が合致する部位はクオリティスコアが上昇し、オーバーラップする領域で塩基が合致しない領域はクオリティスコアが高い方の塩基にbase callが変化したfastqが出力される。マージしたfastq自体は出力されない。クオリティスコアが同じで塩基が合致しない部位はNになる。
・k-merのgraphを利用し、マージされなかったfastqを延長しマージを試みる。
bbmerge-auto.sh in=reads.fq out=merged.fq outu=unmerged.fq ihist=ihist.txt ecct extend2=20 iterations=5
- extend2=0 Extend reads this much only after a failed merge attempt, or in rem/rsem mode.
- iterations=1 (ei) Iteratively attempt to extend by extend2 distance and merge up to this many times.
上記のコマンドだと最初にマージが試みられ、ダメならエラー修復して再チャレンジし、それでもダメならk-merのグラフから20-bpくっつけて再マージを試みる。そうして20-bpのextensionを最大5回繰り返す。途中でk-merのgraphがdead endになったり、分岐があったりする試行は中止される。300-bpx2のペアリードでテストすると、最大794-bpのシングルリードが出力された。面白いコマンドである。
・リード長よりインサートサイズが短いペアリードからアダプター配列を見つける。
bbmerge.sh in1=reads1.fq in2=reads2.fq outa=adapters.fa
例えば見つかったアダプターがGATCGGAAGAGCACACGTCTGAACTCCAGTC とする。このアダプター配列を次のコマンドで除く。
bbmerge.sh in1=reads1.fq in2=read2.fq out=merged.fq adapter1=GATCGGAAGAGCACACGTCTGAACTCCAGTC
リードに十分なデプスがあり、アダプター配列も除かれていれば、mergeコマンドは正しく機能する。より長いリードができればアセンブルでより長いk-merを選択することができることになる。一方で、誤ったgraphをつないでしまう可能性もあり、そうなるとアーティファクトの配列ができてしまうことになる(特にレアなindel変異の部位)。使い方次第で毒にも薬にもなるコマンドかもしれない。
カバレッジがディープな領域をノーマライズして軽量化し、データ解析を行いやすくするツール。例えば、カバレッジが非常に多く変動するシングルセルやウィルスのシーケンスデータ、またメタゲノムデータなどに利用できるとされる。また、シーケンスデータ軽量化のためのdigital normalizationに使えるらしい(リンク 真ん中付近)。リードが多すぎて計算リソースを超えている時などに利用する。
bbnorm.sh in1=reads1.fq in2=reads2.fq out=normalized.fq target=100 min=5
- in=<file> (in1) Input file, or comma-delimited list of files.
- in2=<file> Optional second file for paired reads.
- k=31 Kmer length (values under 32 are most efficient, but arbitrarily high values are supported)
- target=100 (tgt) Target normalization depth. NOTE: All depth parameters control kmer depth, not read depth.
- mindepth=5 (min) Kmers with depth below this number will not be included when calculating the depth of a read.
- maxdepth=-1 (max) Reads will not be downsampled when below this depth, even if they are above the target depth.
- minkmers=15 (mgkpr) Reads must have at least this many kmers over min depth to be retained. Aka 'mingoodkmersperread'.
- out=<file> File for normalized or corrected reads. Use out2 for paired reads in a second file
上記のコマンドではk-merデプスが5以下のリード(シーケンスエラー、またはレアな変異)は排除される。またk-merカバレッジが100になるようにリード数が調整される。
・エラーコレクション。
ecc.sh in1=reads1.fq in2=reads2.fq out=corrected.fq ecc=t keepall passes=1 bits=16 prefilter
- ecc=f Set to true to correct errors.
- keepall=f Set to true to keep all reads (e.g. if you just want error correction).
- passes=2 (p) 1 pass is the basic mode. 2 passes (default) allows greater accuracy, error detection, better contol of output depth.
- bits=32 Bits per cell in bloom filter; must be 2, 4, 8, 16, or 32. Maximum kmer depth recorded is 2^cbits. Automatically reduced to 16 in 2-pass.
- prefilter=f True is slower, but generally more accurate; filters out low-depth kmers from the main hashtable. The prefilter is more memory-efficient because it uses 2-bit cells.
エラーコレクション後のfastqがcorrected.fq として出力される。レアk-merのリードのトリミングや除去は一切行われない。
・k-merヒストグラムの作成。
khist.sh in=reads.fq khist=khist.txt peaks=peaks.txt
ヒストグラムデータの他に、k-mer頻度分布のピークなどをプリントしたサマリーファイルが出力される。
user$ cat peaks.txt
#k 31
#unique_kmers 32597468
#main_peak 14
#genome_size 6316924
#haploid_genome_size 6316924
#fold_coverage 14
#haploid_fold_coverage 14
#ploidy 1
#percent_repeat 36.469
#start center stop max volume
6 14 35 338679 4013227
41 51 121 3023 69382
823 1556 1569 10 1546
1569 1614 1624 21 605
1624 1637 1667 12 455
1667 1683 1705 9 340
1705 1730 1748 13 388
1748 1768 2257 3 1902
3198 4299 4410 4 3095
4410 4501 4523 2 288
4523 4542 5077 2 1112
注意点として、BBnormは発現解析(RNA seqやchip-seq)やレアバリアントの探索などに使ってはならない。 使うとひどいbiasが生じる可能性がある。そのようなデータセットでデータを削減する必要があるなら、ノーマライズでなくランダムサンプリングするツールを使う。
fastqのランダムサンプリングツールとしてseqtkやseqkitを紹介しています。
Dedupe(リンク)-duplicateしたcontigの除去。
k-merベースのde brujin graphでは重複したcontigはほぼ作られないが、overlap layput方式のアセンブルだとduplication contigができることがある。そのようなcontigから重複したcontigを除けるように本コマンドは設計された。また、NCBIのrefseqなどは重複した配列が大量に存在し、blast解析すると同じような配列が大量にhitすることがあるが、そのようなデータベースのキュレーションにも使えるとされる。
また、NCBIの
・duplication sequencesの除去。
dedupe.sh in=X.fa out=Z.fa outd=duplicates.fa
- outd
outdの指定でduplication readsも別出力可能。
・5ミスマッチと2editsまで許容してduplication sequencesを除去。
dedupe.sh in=X.fa out=Y.fa s=5 e=2
s=5だけ指定すると塩基のミスマッチだけ5つまで許容してduplication sequencesが探索される。e=2だけ指定すると、どのような変異でも2つまで許容してduplication sequencesが探索される(例えば2塩基ミスマッチ or 1挿入と1塩基ミスマッチ)。
・99%以上同一の配列を除去。
dedupe.sh in=X.fa out=Y.fa minidentity=99
indelは許容されない。indelも考慮するならe flagを使用する。
・200bp以上同一の配列を1塩基のミスマッチまで許容して除去。
dedupe.sh in=X.fq pattern=cluster%.fq ac=f am=f s=1 mo=200 c pc csf=stats.txt outbest=best.fq fo c mcs=3 cc dot=graph.dot
duplication sequenceというよりアンプリコンシーケンスデータ向けのコマンド。似た配列はクラスタリングして出力され、同時に長さとクオリティがベストの配列は別ファイルに出力される。
・PacBioの16S全長シーケンスデータのクラスタリング。
1、全長読めてない配列を除去。
reformat.sh in=reads_of_insert.fastq out=filtered.fq minlen=1420 maxlen=1640 maq=20 qin=33
- minlength=0 (ml) Reads shorter than this after trimming will be discarded. Pairs will be discarded only if both are shorter.
- maxlength=0 If nonzero, reads longer than this after trimming will be discarded.
- minavgquality=0 (maq) Reads with average quality (after trimming) below this will be discarded.
- qin=auto ASCII offset for input quality. May be 33 (Sanger), 64 (Illumina), or auto.
1420-bp以下のリードと1640-bp以上のリード(アーティファクト)を除去する。また平均クオリティが20以下のリードも除去する。
2、99%以上同一の配列をクラスタリング。
dedupe.sh in=filtered.fq csf=stats_e26.txt outbest=best_e26.fq qin=33 usejni=t am=f ac=f fo c rnc=f mcs=3 k=27 mo=1420 ow cc pto nam=4 e=26 pattern=cluster_%.fq dot=graph.dot
- usejni=f (jni) Do alignments in C code, which is faster, if an edit distance is allowed. This will require compiling the C code; details are in /jni/README.txt.
- absorbmatch=t (am) Absorb exact matches of contigs.
- absorbcontainment=t (ac) Absorb full containments of contigs.
- findoverlap=f (fo) Find overlaps between contigs (containments and non-containments). Necessary for clustering.
- renameclusters=f (rnc) Rename contigs to indicate which cluster they are in.
- ow=f (overwrite) Overwrites files that already exist.
- cc=t (canonicizeclusters) Flip contigs so clusters have a single orientation.
- pto=f (preventtransitiveoverlaps) Do not look for new edges between nodes in the same cluster.
- minclustersize=1 (mcs) Do not output clusters smaller than this.
- cluster=f (c) Group overlapping contigs into clusters.
- numaffixmaps=1 (nam) Number of prefixes/suffixes to index per contig. Higher is more sensitive, if edits are allowed.
- maxedits=0 (e) Allow up to this many edits (subs or indels). Higher is slower.
4塩基オーバーラップの27merのseedを使い、99%以上同一で(e=26)1420塩基以上オーバーラップする配列をクラスタリングする。
コマンドが多いので、2回に分けます。
引用
BBMap (aligner for DNA/RNAseq) is now open-source and available for download.
BBMap (aligner for DNA/RNAseq) is now open-source and available for download. - SEQanswers
BBMerge - Accurate paired shotgun read merging via overlap.
Bushnell B, Rood J, Singer E
PLoS One. 2017 Oct 26;12(10):e0185056. doi: 10.1371/journal.pone.0185056. eCollection 2017.
BIostars ダウンサンプリング
https://www.biostars.org/p/168178/