多機能なNGSの管理ツール BBtools　其の１ - macでインフォマティクス

2018 9/5 bbmerge-auto.sh修正

BBtoolsはアメリカのJGIが提供している多機能なNGS向けの解析ツール。2014年にオープンソース化されたらしい。論文は現在準備中とある。アライメントのBBmapや、オーバーラップがないペアリードをマージするBBMerge、エラーコレクションしたfastqを出力するBBnormなど特徴的なコマンドを備える。ここでは以下のコマンドを紹介する。

BBMap - リファンレンスへのアライメント。
BBMerge - ペアリードのマージ。k-mer配列を付け足し延長することも可能。
BBnorm - カバレッジのノーマライズ。
Dedupe - duplicateしたcontigの除去。

其の２。

公式ページ

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 （リンク）-リードのリファンレンスへのアライメント

powerpointスライド

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.

さらに指定した分析データ（ヒストグラム）が出力される。

f:id:kazumaxneo:20170829140136j:plain

参考１

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の物理メモリが必要と算出された。

参考２

indel部位のスプリットアライメントに対応している。条件はmaxindelフラッグで規定されており、デフォルトではmaxindel=16000となっている。RNA seqでlarge intronが想定される場合などはmaxindel=100000などの大きな数値に指定することが推奨されている。

BBMerge （リンク）-ペアリードのマージ。

BBMegeはオーバーラップするペエアリードをマージして１つの長いリードを作るツール。

・ペアリードのマージ。

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 in1=reads1.fq in2=reads2.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を最大５回繰り返す。途中で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変異の部位）。使い方次第で毒にも薬にもなるコマンドかもしれない。

BBnorm （リンク）-カバレッジのノーマライズ。

カバレッジがディープな領域をノーマライズして軽量化し、データ解析を行いやすくするツール。例えば、カバレッジが非常に多く変動するシングルセルやウィルスのシーケンスデータ、またメタゲノムデータなどに利用できるとされる。また、シーケンスデータ軽量化のためのdigital normalizationに使えるらしい（リンク　真ん中付近）。リードが多すぎて計算リソースを超えている時などに利用する。

・k-merカバレッジ100でノーマライズ。

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だけ指定すると塩基のミスマッチだけ５つまで許容してduplication sequencesが探索される。e=2だけ指定すると、どのような変異でも２つまで許容してduplication sequencesが探索される（例えば2塩基ミスマッチ or 1挿入と1塩基ミスマッチ）。

・99%以上同一の配列を除去。

dedupe.sh in=X.fa out=Y.fa minidentity=99

indelは許容されない。indelも考慮するならe flagを使用する。

・200bp以上同一の配列を１塩基のミスマッチまで許容して除去。

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.