2018 11/22 ポスターlink追加、誤字修正
次世代シーケンシング技術（NGS）は、シーケンシングされたゲノムの生成を大幅に加速した。しかしながら、これらの技術は、依然として染色体当たり単一の配列を提供することができないままである。代わりに、それらは大量かつ冗長なリードセットを生成し、各リードは全ゲノムの一部である。この冗長性のために、リード間のオーバーラップを検出し、それらを一緒にアセンブルしてターゲットのゲノム配列を再構成することが可能である。
今日でさえ、リードのアセンブリは複雑な作業のままであり、これはどのようなソフトウェアも一貫してうまく働かないためである[論文より ref.1]。アセンブリの問題は以前からずっと計算上困難なままであり、より正確にはNP困難であることが示されている[ref.2]。実用上の限界は、ゲノムの構造（正確に解読できるより長いリピート）およびシーケンシングバイアス（不均一なカバレッジおよびシークエンシングエラー）の両方から生じている。適用されたソリューションは、リードをアセンブリgraphとして表現する: graphのパスに沿ったラベルがシーケンスをエンコードする。現在のところ、大部分のアセンブラは、2種類のgraphに依存している。第2世代のシーケンシング技術によって生成されたショートリードのde Bruijn graph（DBG）[ref.3]、またはCelera Assembler [ref.4]）とそのバリエーション、string graph [ref.5]のように。次に、アセンブリアルゴリズムはヒューリスティックを使用してgraphを探索し、いくつかの経路を選択してそのシーケンスを出力する。これらのヒューリスティックのために、得られた配列の集合（コンティグと呼ばれる）は、シーケンシングエラー、ゲノムの変異およびリピートによって、生成されるgraphは複雑なパターンになり断片化される。コンティグセットはめったに満足のいくものではなく、通常短いコンティグを破棄するなどして後処理される。
リードセットを分析するための最も頻繁なコンピュータタスクは、それらをリファレンスゲノム上にマッピングすることである。リファレンスゲノムが一連の配列の形の場合、リードをマップすることができるツールはたくさんある（例えば、BWA [ref.6]およびBowtie [ref.7]）。Finishedゲノム配列にマッピングすることの目的は、アライメントできるかどうか、そしてできる場合はどのポジションにアライメントされるか、である。これには、たいてい、リード配列とゲノム配列の間で小さな編集距離またはハミング距離を許可するsemi-globalなアラインメント手順のヒューリスティックで行われる。リードマッピングプロセスは、低マッピング率の領域に苦しんでいる[ref.8]。リピートゲノム領域では、リードのマッピングが複数一致するため、正確にマッピングされないかもしれない。ゲノムがgraphとして表されるとき、各リピート領域出現は分解されて複数一致問題は制限されるため、マッピング問題は低減される。
リファレンスがFinishedゲノム配列ではなく、重複したコンティグセットである場合、状況は異なる。マッピングは、リードがゲノム内に見出されるかどうかを正確に決定し得るが、複数の位置は、例えば、いくつかの真の位置が存在するかどうかを判断するには不十分であり得る。逆に、不完全なマッピングは、不完全なアセンブリまたは後処理中のいくつかのコンティグの除去に起因する可能性がある。そのような場合、コンティグセットの代わりに、アセンブリgraphを（偏っていないおよび/またはより完全な）リファレンスとして考えることは興味深いと我々（著者ら）は主張する。次いで、このグラフ経路のマッピングは、コンティグのセットに対するマッピングを補完するために必要である。これがBGREATの設計と実装を促した。
ここでは、graph上にリードをマッピングする問題を探究する。リードをsequence graph（ (a generic term meaning a graph representing sequences along its paths) ）にマッピングまたはアライメントすることは、既にいくつかの文献、例えばアセンブリ、リードのエラー訂正、メタゲノミクス分野で検討されてきたことである。
アセンブリ関連では、一旦DBGが構築されると、リードをgraphにバックマッピングすることで、サポートされていないパスを削除したり、エッジのカバレッジを計算したりすることができる。著者らが知る限り、この作業のための実用的な解決法は設計されていない（論文執筆時点）。 Ceruleanアセンブラ[ref.9]はこの可能性について言及しているが、通常のアセンブリ配列へのアラインメントしか実行しない。 Allpaths-LG [ref.10]も、第3世代シーケンシング技術のエラーの多いロングリードを使用してリピートを解決する同様のタスクを実行するが、その手順は、DBG上の任意のリードセットのマッピングに適するほど一般的ではない。理論的な観点からは、問題は、NP困難なリード通過問題（オイラーのスーパーパス問題[2,11]とも呼ばれる）であり、DBGでリードのコヒーレントパスを見つけることからなる（ [ref.5]で定義されているような一連のリード配列として表現できるパス）。 SPADES [ref.12]アセンブラは、構築中に使用されるパスを追跡することによってDBGに対するリードを通す。これには相当量のメモリが必要である。ここでは、グラフの作成に使用されるリードやその他のリードのいずれかをNGSリードの任意のソースにマップすることを目指して、より一般的な問題を提案する（De Bruijn Graph Read Mapping Problem（DBGRMPと名付ける））。
最近では、ショートリードを使ったロングリードのハイブリッドエラー訂正は、第3世代のシーケンシング技術を活用するための重要なステップとなっている。エラー訂正ツールLoRDEC [ref.13]は、ショートリードのDBGを構築し、動的プログラミングアルゴリズム（著者らの目的にとっては遅い）を用いてそれらの編集距離を計算することによってDBGの各行に対して各ロングリードをアライメントする。ショートリードの訂正ツールでは、シーケンシングエラーを訂正するために評価するためにk-merスペクトラムを使うツールが、DBGの確率的または正確な表現をリファレンスとして使ってる[ref.14,15]。
メタゲノミクス関連では、Wang et al [ref.16]が、closely relatedな複数の細菌種からなるDBG上にリードをマッピングすることによって、サンプルの分類組成を推定した。実際、graphはこれらのゲノムの類似領域を崩壊させ、冗長マッピングを回避する。彼らのツールは、DBGのunitigsのランダムな連結配列にBWAを使用してリードをマッピングする。したがって、graphのいくつかの連続したノード（ER: il y a un pb ce n’est pas vrai)）上にリードをアライメントすることはできない。同様に、いくつかの著者は、関連するゲノムを、より反復性の少ない単一のDBGに保存することを提案している[17,18,19]。しかし、これらのツールのほとんどは、hihgly closely relatedな配列に適用した場合にのみ効率的であり、結果としてフラットなgraphになる。 BlastGraphツール[19]は、graphのリードマッピングに特化しているが、実際のgraphでは使用できない（論文の結果セクションを参照）。
ここでは、De Bruijn graph上のリードマッピングを正式化し、それがNP完全であることを示す。次に、パイプラインGGMAPを提示し、DBGの分岐パス上のリードをマップできる新しいツールBGREAT (Section GGMAP: a method to map reads on de Bruijn Graph)に入る。効率性のために、BGREATは、巨大なシーケンシングデータセットまでスケールアップするヒューリステイックスを採用している。セクションの結果では、マッピング能力と効率の面でGGMAPを評価し、それをアセンブルされたコンティグのマッピングと比較する。最後に、GGMAPの限界と利点について議論し、今後の作業の方向性について述べる（Section Discussion）。
ポスター
インストール
ubuntu14.04でテストした（docker使用。ホストOS: mac os10.12）。
本体 Github
git clone https://github.com/Malfoy/BGREAT.git
cd BGREAT/
make
>./bgreat
$ ./bgreat
-r read_file
-k k_value (30)
-g unitig_file (unitig.dot)
-m n_missmatch (2)
-t n_thread (1)
-e effort put in mapping (2)
-f path_file (paths)
-a not_aligned_file (notAligned.fa)
-p to align on paths instead of walks
-q for fastq read file
-c to output corrected reads
BCALM (constructs the compacted de Bruijn graph from sequencing data)
#v1 (DSKも必要になる）
git clone https://github.com/Malfoy/BCALM.git
cd BCALM/
make
#またはcondaで導入(version2.2)(github)
conda install -y -c bioconda bcalm
>./bcalm #version1
$ ./bcalm
usage: <input> [output.dot] [minimizer length]
Note: maximum minimizer length (minimizer length = 10) requires that you type 'ulimit -n 1100' in your shell prior to running bcalm, else the software will crash
> bcalm #version2
$ bcalm
BCALM 2, git commit c8ac60252fa0b2abf511f7363cff7c4342dac2ee
Specifiy -in
[bcalm_1 options]
-nb-cores (1 arg) : number of cores [default '0']
-verbose (1 arg) : verbosity level [default '1']
-version (0 arg) : version
-help (0 arg) : help
[graph options]
-no-mphf (0 arg) : don't construct the MPHF
[kmer count options]
-in (1 arg) : reads file [default '']
-kmer-size (1 arg) : size of a kmer [default '31']
-abundance-min (1 arg) : min abundance threshold for solid kmers [default '2']
-abundance-max (1 arg) : max abundance threshold for solid kmers [default '2147483647']
-solidity-custom (1 arg) : when solidity-kind is cutom, specifies list of files where kmer must be present [default '']
-max-memory (1 arg) : max memory (in MBytes) [default '5000']
-max-disk (1 arg) : max disk (in MBytes) [default '0']
-out (1 arg) : output file [default '']
-out-dir (1 arg) : output directory [default '.']
-out-tmp (1 arg) : output directory for temporary files [default '.']
-out-compress (1 arg) : h5 compression level (0:none, 9:best) [default '0']
-storage-type (1 arg) : storage type of kmer counts ('hdf5' or 'file') [default 'hdf5']
[kmer count, algorithmic options options]
-minimizer-type (1 arg) : minimizer type (0=lexi, 1=freq) [default '1']
-minimizer-size (1 arg) : size of a minimizer [default '10']
-repartition-type (1 arg) : minimizer repartition (0=unordered, 1=ordered) [default '1']
[bloom options]
-bloom (1 arg) : bloom type ('basic', 'cache', 'neighbor') [default 'neighbor']
-debloom (1 arg) : debloom type ('none', 'original' or 'cascading') [default 'cascading']
-debloom-impl (1 arg) : debloom impl ('basic', 'minimizer') [default 'minimizer']
[branching options]
-branching-nodes (1 arg) : branching type ('none' or 'stored') [default 'stored']
-topology-stats (1 arg) : topological information level (0 for none) [default '0']
[general options]
-config-only (0 arg) : dump config only
-nb-cores (1 arg) : number of cores [default '0']
-verbose (1 arg) : verbosity level [default '1']
-integer-precision (1 arg) : integers precision (0 for optimized value) [default '0']
[debug options]
-redo-bcalm (0 arg) : debug function, redo the bcalm algo
-skip-bcalm (0 arg) : same, but skip bcalm
-redo-bglue (0 arg) : same, but redo bglue
-skip-bglue (0 arg) : same, but skip bglue
-redo-links (0 arg) : same, but redo links
-skip-links (0 arg) : same, but skip links
DBGを出力するならDSKなどのk-merカウンターも必要になる。
conda install -y -c bioconda dsk
実行方法
1、k-merカウント及びde Bruijn graph構築（BCALM v1）
#BCALM v1を使う場合
#k-merカウント
dsk -file test.fa -kmer-size 31 -abundance-min-threshold 3 -out out
dsk2ascii -file out.h5 -out k-mer_count
bcalm k-mer_count output.dot 10
> k-mer_counthead
$ k-mer_counthead
bash: k-mer_counthead: command not found
kazu@6a393ab918c4:/data$ head k-mer_count
AAAAAAAGATTTGATTGCGCCGATGGTTTTT 23
AAAAAATCATTGGCACGTTTGGCAACTGTGC 9
AAAAAAGATTTGATTGCGCCGATGGTTTTTA 23
AAAAAAGTGCGCACATGTGCGCACTTTTCTG 12
AAAAACAGCCTGAACTATATGACTTTAACGG 15
AAAAACCAGCGGCGATCGTCCTTAACCTACT 7
AAAAACGTGGAGATCCTCAGGAGTTTTACCT 14
AAAAATCATTGGCACGTTTGGCAACTGTGCT 9
AAAAAGATTTGATTGCGCCGATGGTTTTTAC 23
AAAAAGTGCGCACATGTGCGCACTTTTCTGC 12
or
de Bruijn graph構築（BCALM v2 Github)
bcalm -in input.fq -kmer-size 21 -abundance-min 5
#複数のシーケンスデータがあるならリストで与える
ls -1 *.fastq > list_reads
bcalm -in list_reads -kmer-size 21 -abundance-min 10
version2では、unitigのfastaファイルが出力される(*1 GFA変換)。
２、 de Bruijn graphにマッピング。シーケンスデータをfasta形式、またはfastq形式（"-q"をつける）で指定する。
bgreat -r reads.fa,reads2.fa -k 21 -t 12 -m 1 -g unitig.fa \
-o no_overlap -a not_aligned -p path
出力の詳細はGithubを参照してください。
引用
Read mapping on de Bruijn graphs
Limasset A, Cazaux B, Rivals E, Peterlongo P
BMC Bioinformatics. 2016 Jun 16;17(1):237
*1
付属のツールを使い、GFAフォーマット（下にリンク）に変換できる。21はDBG構築時に指定したk-merサイズ。
python2 bcalm/scripts/convertToGFA.py -s unitigs.fa output.gfa 21
出力
GFAファイルはbandageに読み込んで可視化できる。
The GFA Format Specification
https://github.com/GFA-spec/GFA-spec/blob/master/GFA1.md
