20196/13 タイトル修正
2020 7/24 bbdukコマンド追記
の続き。BBtoolsの残りのコマンドを紹介する。紹介するのは以下のコマンド。
- Reformat - フォーマット変換やクオリティトリミング。
- Repair - ペアリードの順番が壊れたファイルを修復する。
- Stats - アセンブリの基本情報をレポートする。
- BBDuk - クオリティやアダプターでトリミングしたり、フィルタリングする。
- BBMask - 単純リピートや保存性が高い領域をマスクする。
-
CalcUniqueness - ユニークなリードの数を調べる。
インストール
https://sourceforge.net/projects/bbmap/
#bioconda (link)
mamba install -c bioconda -y bbmap
#homebrew
brew install BBtools
Reformat（リンク）-フォーマット変換やクオリティトリミング。
似たコマンドはいくつかあるが、Reformatはメモリー使用量が少ないのが特徴。
下では全てインターレースのfastqを想定しているが、ペアリードなら
reformat.sh in=<file> in2=<file2> out=<outfile> out2=<outfile2> と書く。
・fastq => fastaへの変換。
reformat.sh in=reads1.fq out=reads1.fasta
・phread quality scoreを一般的なASCII+33からASCII +64に変換する。
reformat.sh in=reads.fq out=reads.fq qin=33 qout=64
・fastq => fasta+qualへの変換。
reformat.sh in=reads.fq out=reads.fa qfout=reads.qual
・ インターレースのペアリードデータを２つのファイルに分ける。
reformat.sh in=interlace.fq out1=read1.fq out2=read2.fq
・ ペアリードデータをインターレースファイルにする。
reformat.sh in1=read1.fq in2=read2.fq out=reads.fq
・ 70文字ごとに改行する。(foldコマンドと同じ）
reformat.sh in=reads.fa out=wrapped.fa fastawrap=70
・ペアリードファイルがずれていないか検証。
reformat.sh in=reads.fq vint
- verifyinterleaved=f (vint) sets 'vpair' to true and 'interleaved' to true.
・1や２がついていないペアリードファイルにペアを識別する1 or 2を入れる。
reformat.sh i1n=pair1.fq in2=pair2.fq out1=renamed1.fq out1=renamed2.fq addslash=t int=f
- addslash=f Append ' /1' and ' /2' to read names, if not already present. Please include the flag 'int=t' if the reads are interleaved.
- int=f (interleaved) Determines whether INPUT file is considered interleaved.
・reverse complement
reformat.sh in=reads.fq out=out.fq rcomp
- rcomp=f (rc) Reverse-compliment reads.
・小文字の塩基を大文字にする。
reformat.sh in=reads.fq out=out.fq tuc
- touppercase=t (tuc) Convert input bases to upper-case; otherwise, lower-case will not match.
・名前が全てユニークか確認する。
reformat.sh in=reads.fq out=out.fq uniquenames
- uniquenames=f Make duplicate names unique by appending _<number>.
duplicationは"_number"がついてユニークな名前になる。このコマンドだけ大量のメモリを使う。
・クオリテイスコアを特定の範囲に限定する。
reformat.sh in=reads.fq out=out.fq mincalledquality=2 maxcalledquality=41
- maxcalledquality=41 Quality scores capped at this upper bound.
- mincalledquality=2 Quality scores of ACGT bases will be capped at lower bound.
特定のツールの処理でクオリティがおかしくなってしまった時などに使う。
Repair（リンク）-ペアリードの順番が壊れたファイルを修復する。
ペアリードの順番がforwardとreverseでずれているfastqを修復する。例えばfastx toolkitでペアリードを個別にトリミングすると、forwardとreverseでリード数が揃わないfastqができてしまう。この場合、理想はraw readsからペアリードを扱えるツールでクオリティトリミングをやり直すことであるが、場合によってはオリジナルのfastqがないこともある。そのようなfastqの修復にこのコマンドを使う。
入力 - fasta、fastq、sam。
・ペアリードファイルの修復。
repair.sh -Xmx20g in1=broken1.fq in2=broken2 out1=fixed1.fq out2=fixed2.fq outs=singletons.fq repair
- in=<file> The 'in=' flag is needed if the input file is not the first parameter. 'in=stdin' will pipe from standard in.
- in2=<file> Use this if 2nd read of pairs are in a different file.
- out=<file> The 'out=' flag is needed if the output file is not the second parameter. 'out=stdout' will pipe to standard out.
- out2=<file> Use this to write 2nd read of pairs to a different file.
- outs=<file> (outsingle) Write singleton reads here.
- -Xmx20g This will be passed to Java to set memory usage, overriding the program's automatic memory detection. -Xmx20g will specify 20 gigs of RAM, and -Xmx200m will specify 200 megs. The max is typically 85% of physical memory.
・contig.fasstaを分析する。
stats.sh in=contigs.fa gc=gc.txt gcformat=4
- in=file Specify the input fasta file, or stdin.
- gc=file Writes ACGTN content per scaffold to a file.
- gcformat=<0-4> Select GC output format; default 1.
- gcformat=0: (no base content info printed)
- gcformat=1: name length A C G T N GC
- gcformat=2: name GC
- gcformat=4: name length GC
Note that in gcformat 1, A+C+G+T=1 even when N is nonzero.
N50などの基礎的な情報が出力される。
user$ stats.sh in=contigs.fa gc=gc.txt gcformat=4
A C G T N IUPAC Other GC GC_stdev
0.2841 0.2276 0.2531 0.2352 0.0000 0.0000 0.0000 0.4807 0.0583
Main genome scaffold total: 42
Main genome contig total: 42
Main genome scaffold sequence total: 0.004 MB
Main genome contig sequence total: 0.004 MB 0.000% gap
Main genome scaffold N/L50: 9/131
Main genome contig N/L50: 9/131
Main genome scaffold N/L90: 42/41
Main genome contig N/L90: 42/41
Max scaffold length: 352
Max contig length: 352
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 42 42 3,520 3,520 100.00%
50 19 19 2,547 2,547 100.00%
100 11 11 2,009 2,009 100.00%
250 1 1 352 352 100.00%
gc含量は別ファイルに出力される。
user$ head -5 gc.txt
#Name Length GC
NODE_1_length_332_cov_90.186745 352 0.5114
NODE_4_length_46_cov_61.913044 66 0.4242
NODE_6_length_127_cov_8.811024 147 0.3469
NODE_7_length_127_cov_14.677165 147 0.2993
BBDuk（リンク）-クオリティトリミング、アダプタートリミング、フィルタリング、マスキングなどを行う。
BBDukはクオリティトリミングを行ったり、k-merベースでコンタミリードと推測されるリードのフィルタリングを行うことができるツール。GC、サイズ、エントロピーなど多様な方法でフィルタリングすることもできる。
ペアリードの場合、
bbduk.sh in1=read1.fq in2=read2.fq out1=clean1.fq out2=clean2.fq
というような書き方をする。BBDukは常にペアリードの順番が壊れないように動作する。正常動作のためにペアリードは必ず同時に処理する。
・3'側のアダプタートリミング。
bbduk.sh in=reads.fq out=clean.fq ref=adapters.fa ktrim=r k=23 mink=11 hdist=1 tpe tbo
- in=<file> Main input. in=stdin.fq will pipe from stdin.
- out=<file> (outnonmatch) Write reads here that do not contain kmers matching the database. 'out=stdout.fq' will pipe to standard out.
- ref=<file,file> Comma-delimited list of reference files.
- tpe=f (trimpairsevenly) When kmer right-trimming, trim both reads to the minimum length of either.
- hammingdistance=0 (hdist) Maximum Hamming distance for ref kmers (subs only). Memory use is proportional to (3*K)^hdist.
- tbo=f (trimbyoverlap) Trim adapters based on where paired reads overlap.
- mink=0 Look for shorter kmers at read tips down to this length, when k-trimming or masking. 0 means disabled. Enabling this will disable maskmiddle.
- ktrim=f Trim reads to remove bases matching reference kmers. Values:
f (don't trim), f (don't trim),
r (trim to the right),
l (trim to the left)
ktrim=rでright trimming（3'側アダプター）、ktrim=lでleft trimming（5'側アダプター）を行う。kでアダプターのを探すサイズ（アダプターサイズにする）、minkで末端でのサイズを規定する。hdistで
・3'側のクオリティが10以下の領域をトリミング
bbduk.sh in=reads.fq out=clean.fq qtrim=r trimq=10
- trimq=6 Regions with average quality BELOW this will be trimmed.
- qtrim=f Trim read ends to remove bases with quality below trimq. Performed AFTER looking for kmers. Values:
rl (trim both ends),
f (neither end),
r (right end only),
l (left end only),
w (sliding window).
qtrim=rで3'側のクオリティトリミングが実行される。qtrim=rlで両側トリミング。qtrim=f以外のflagが指定されるとトリミングが実行されるが、トリミングはk-merベースの作業が全て完了してから実行される。
・5'側10塩基と3'側139塩基以降を強制トリミング。
bbduk.sh in=reads.fq out=clean.fq ftl=10 ftr=139
- forcetrimleft=0 (ftl) If positive, trim bases to the left of this position (exclusive, 0-based).
- forcetrimright=0 (ftr) If positive, trim bases to the right of this position (exclusive, 0-based).
・3'末端を強制トリミング。
bbduk.sh in=reads.fq out=clean.fq ftm=5
- forcetrimmod=0 (ftm) If positive, right-trim length to be equal to zero, modulo this number.
５で剰余演算（モジュロ）を行い、余りの1がトリムされる。イルミナのシーエンスは基本5単位で増減し、50、75、150、300などがあるが、全て1-bp余分に読まれる。このラストの1-bpのクオリティが非常に低い傾向にあり、そのため5で割って余りの値1でトリムしている。
・平均クオリティが10以下のリードをフィルタリング（除去）。
bbduk.sh in=reads.fq out=clean.fq maq=10
- minavgquality=0 (maq) Reads with average quality (after trimming) below this will be discarded.
・平均クオリティ10以下の領域を トリムし、100bp以下になったリードは除去する。
bbduk.sh in=reads.fq out=clean.fq qtrim=r trimq=10 ml=100
- trimq=6 Regions with average quality BELOW this will be trimmed.
- minlength=10 (ml) Reads shorter than this after trimming will be discarded. Pairs will be discarded if both are shorter.
・k-merフィルタリング。
bbduk.sh in=reads.fq out=unmatched.fq outm=matched.fq ref=phix.fa k=31 hdist=1 stats=stats.txt
#paied-end ( マップもアンマップもペアは同期される)
bbduk.sh in1=pair_1.fq in2=pair_2.fq \
out1=unmatched_1.fq out2=unmatched_2.fq \
outm1=matched_1.fq outm2=matched_2.fq \
ref=phix.fa k=31 hdist=1 stats=stats.txt
- outm=<file> (outmatch) Write reads here that contain kmers matching the database.
- ref=<file,file> Comma-delimited list of reference files.
- hammingdistance=0 (hdist) Maximum Hamming distance for ref kmers (subs only). Memory use is proportional to (3*K)^hdist.
指定ファイルに対して31-merで探索を行い、1塩基のミスマッチを許容して（hdist=1）マッチするリードがあれば、除去される。上ではイルミナのスパイク（バクテリオファージ、5.5-kb）に対して検索をかけている。ペアエンドはマッチしたものもマッチしなかったものも同期される。
動植物ゲノムアセンブリ前に微生物由来リードを除く。
・k-merマスキング１。
bbduk.sh in=ecoli.fa out=ecoli_masked.fa ref=salmonella.fa k=25 ktrim=N
25-merでサルモネラゲノムとマッチするリードを探し、マッチする部位がNNNでマスクされる。Xでマスクするならktirm=X。
・k-merマスキング２。エントロピーを指定。
bbduk.sh in=reads.fq out=filtered.fq entropy=0.5 entropywindow=50 entropyk=5
- entropy=-1 Set between 0 and 1 to filter reads with entropy below that value. Higher is more stringent.
- entropywindow=50 Calculate entropy using a sliding window of this length.
- entropyk=5 Calculate entropy using kmers of this length.
平均エントロピーが0.5以下のリードを除去する。エントロピーは塩基の複雑さから判定している。AAAAAAAAAAAAAAのような配列のエントロピーは０で、ランダムな配列ほど1に近くなる。よって単純リピートなども排除されると考えられる。
・Quality recalibration
bbduk.sh in=reads.fq out=recalibrated.fq recalibrate sam=mapped.sam
- recalibrate=f (recal) Recalibrate quality scores. Requires calibration matrices generated by CalcTrueQuality.
- sam=<file,file> If recalibration is desired, and matrices have not already been generated, BBDuk will create them from the sam file.
samファイルを元にベースクオリティスコアを再校正したfastqを出力する。samのcigar stringは=（マッチ）とX（ミスマッチ）でないといけない。Mとかだと機能しないので、BBMapなどでアライメントをやり直す。リファレンスがない場合、アセンブルしたcontigにアライメントする方法が推奨されている。
・Histogram generation
bbduk.sh in=reads.fq bhist=bhist.txt qhist=qhist.txt gchist=gchist.txt aqhist=aqhist.txt lhist=lhist.txt gcbins=auto
- gcbins=100 Number gchist bins. Set to 'auto' to use read length.
base frequency、クオリティスコア、GC、 平均クオリティスコア、長さ、などの統計情報のヒストグラムデータを出力する。BBMapにも同様のコマンドが存在する。
・Barcode and chastity filtering
bbduk.sh in=reads.fq out=filtered.fq barcodes=ACGTT,GGCCG barcodefilter chastityfilter
- barcodes= Comma-delimited list of barcodes or files of barcodes.
Header-parsing parameters - these require Illumina headers:
- chastityfilter=f (cf) Discard reads with id containing ' 1:Y:' or ' 2:Y:'.
- barcodefilter=f Remove reads with unexpected barcodes if barcodes is set, or barcodes containing 'N' otherwise. A barcode must be the last part of the read header.
・Cardinality estimation
bbduk.sh in=reads.fq loglog loglogk=31
- cardinality=f (loglog) Count unique kmers using the LogLog algorithm.
- loglogk=31 Use this kmer length for counting.
ユニークな31-merの配列の数をカウント。
・Matching degenerate sequences such as primers
bbduk.sh in=reads.fq out=matching.fq literal=ACGTTNNNNNGTC copyundefined k=13 mm=f
- literal=<seq,seq> Comma-delimited list of literal reference sequences.
- copyundefined=f (cu) Process non-AGCT IUPAC reference bases by making all possible unambiguous copies. Intended for short motifs or adapter barcodes, as time/memory use is exponential.
degenerateな配列を探索する。例えばメタ16Sの可変性領域を網羅的に探すことなどに使うことができる。
BBMask（リンク）-単純リピートや保存性が高い領域をマスクする。
false positiveのマッピングを防止するため、単純リピートなどをあらかじめマスクしておくなどの目的に使われる。リードが単純な配列なのか判別するためにシャノン情報量の概念が利用されているらしい。マスクでなくフィルタリングで除くにはBBDukを使用する。
・エントロピー（情報量）が0.7以下の部位をマスクする。
bbmask.sh in=ref.fa out=masked.fa entropy=0.7
- entropy=0.70 (e) Mask windows with entropy under this value (0-1). 0.0001 will mask only homopolymers and 1 will mask everything.
閾値をどのような値にするか正確に求めることは難しいが、entropy=0.7ではE.coliゲノムではわずか70bpしかマスクされない。またコピーや繰り返しが多いヒューマンでも0.7%しかマスクされないとされる。 entropy=0.9まであげるとE.coliで8-kb、ヒューマンで7%の領域がマスクされる。
CalcUniqueness（リンク）-ユニークなリードの数を調べる。
リードの先頭からk-merの配列を取り出しユニークなのかどうか判別する。k-merの配列が合致するとユニークではないと判断される。データ中のユニークなリード数の割合を求めることで、シーケンスのカバレッジが十分かどうか推定することができる。
・ユニークなリード数を示したヒストグラムをかく。
bbcountunique.sh in=reads out=histogram.txt
25000ごとに区切ってユニークなリードの割合が出力される。
k-merの配列は完全に一致していないと重複とみなされない。そのため低クオリティのデータに使うと、見かけ上ユニークなリードの割合が増える点に注意が必要。
他にもいくつかコマンドがあります。公式サイトで確認してみてください。
BBTools User Guide - DOE Joint Genome Institute
引用
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
