macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

多機能なNGS分析ツール BBtools 其の2

20196/13 タイトル修正

2020 7/24  bbdukコマンド追記

 

の続き。BBtoolsの残りのコマンドを紹介する。紹介するのは以下のコマンド。

  • Reformat フォーマット変換やクオリティトリミング。
  • Repair - ペアリードの順番が壊れたファイルを修復する
  • Stats アセンブリの基本情報をレポートする
  • BBDuk - クオリティやアダプターでトリミングしたり、フィルタリングする。
  • BBMask - 単純リピートや保存性が高い領域をマスクする 
  • CalcUniqueness - ユニークなリードの数を調べる。

     

 

 

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

 

インターレースのペアリードデータを2つのファイルに分ける。

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や2がついていないペアリードファイルにペアを識別する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.

 

 

 

Statsリンク)-アセンブリの基本情報のレポート

・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.

 5で剰余演算(モジュロ)を行い、余りの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マスキング1。

bbduk.sh in=ecoli.fa out=ecoli_masked.fa ref=salmonella.fa k=25 ktrim=N

25-merでサルモネラゲノムとマッチするリードを探し、マッチする部位がNNNでマスクされる。Xでマスクするならktirm=X。

 

・k-merマスキング2。エントロピーを指定。

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のような配列のエントロピーは0で、ランダムな配列ほど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