シーケンシング実験あたりのデータ量が増加するにつれて、フラグメントアセンブリはますます計算量が増加している。De Bruijn graphは、フラグメントアセンブリアルゴリズムで広く使用されているデータ構造で、リードのセットからの情報を表現するために使用されている。Compactionは、長い単純なパスを単一の頂点に圧縮するde Bruijn graphベースのアルゴリズムにおいて、データ削減の重要なステップである。最近ではアセンブリパイプラインのボトルネックとなっている圧縮処理の実行時間とメモリ使用量を改善することが重要な課題となっている。
ここではDe Bruijn graphの圧縮のためのアルゴリズムとツールBCALM 2を紹介する。BCALM2はミニマムハッシュ法に基づいて入力を分散させる並列アルゴリズムであり、実行中のメモリ使用量のバランスを考慮している。ヒトのシーケンシングデータに対しては、BCALM2を適用することで、de Bruijn graphの圧縮にかかる計算負荷を約1時間、メモリ使用量を3GBに削減することができた。また、22 Gbpのloblolly pineと20Gbpのwhite spruceのシーケンシングデータにもBCALM 2を適用した。コンパクト化されたグラフは、1台のマシン上で2日以内に生のリードから40GBのメモリで構築された。そのため、BCALM 2は他の方法と比較して、少なくとも一桁以上の効率性がある。
インストール
macos10.14でビルドしてテストした。
ビルド依存
#bioconda (link)
conda install -c bioconda bcalm -y
#from source
git clone --recursive https://github.com/GATB/bcalm
cd bcalm
mkdir build; cd build
cmake ..; make -j 8
#large k-merを扱うなら
git clone --recursive https://github.com/GATB/bcalm
cd bcalm
mkdir build; cd build
rm -Rf CMake* && cmake -DKSIZE_LIST="32 64 96 128 160 192 224 256 320" .. && make -j 8
> ./bcalm
$ ./bcalm
BCALM 2, version v2.2.3, git commit 9b7b581
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 custom, 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']
-histo2D (1 arg) : compute the 2D histogram (with first file = genome, remaining files = reads) [default '0']
-histo (1 arg) : output the kmer abundance histogram [default '0']
[kmer count, advanced performance tweaks 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']
-all-abundance-counts (0 arg) : output all k-mer abundance counts instead of mean
-edge-km (1 arg) : edge km representation [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 (needs debug_keep_glue_files=true in source code)
-skip-bglue (0 arg) : same, but skip bglue
-redo-links (0 arg) : same, but redo links
-skip-links (0 arg) : same, but skip links
-nb-glue-partitions (1 arg) : number of glue partitions
実行方法
fastqのリストファイルを指定する。
ls -1 *.fastq > list_reads
bcalm -in list_reads -kmer-size 55 -abundance-min 2 -out-dir outdir
- -kmer-size size of a kmer [default '31']
- -abundance-min min abundance threshold for solid kmers [default '2']
- -out-dir output directory [default '.']
unitigが出力される。
出力について
https://github.com/GATB/bcalm/blob/master/bidirected-graphs-in-bcalm2/bidirected-graphs-in-bcalm2.md
付属のスクリプトでunitigをGFAに変換する。unitig、出力GFA、k-mer値の順番に指定する。
python bcalm/scripts/convertToGFA.py list_reads.unitigs.fa output.gfa 55
引用
Compacting de Bruijn graphs from sequencing data quickly and in low memory
Rayan Chikhi, Antoine Limasset, Paul Medvedev
Bioinformatics, Volume 32, Issue 12, 15 June 2016, Pages i201–i208
