vg日本語版ドキュメント # はじめに vgの日本語版ドキュメントをまとめると便利だという話になったので、vgに関する情報をまとめていく。とりあえず最近使用したコマンドを記載する。詳しい説明やその他のコマンドについては順次追加していく。 # vgのビルド ## 事前に必要なライブラリ ``` sudo apt-get install build-essential git cmake pkg-config libncurses-dev libbz2-dev \ protobuf-compiler libprotoc-dev libjansson-dev automake libtool \ jq bc rs curl unzip redland-utils librdf-dev bison flex lzma-dev \ liblzma-dev liblz4-dev ``` リポジトリの取得 ``` git clone --recursive https://github.com/vgteam/vg.git ``` ``` . ./source_me.sh && make static ``` でビルド ``` ./bin/vg ``` で実行 ## Mac OS Xでのビルド # vgの使い方 ## vgのコマンドラインインターフェース ``` construct: graph construction view: conversion (dot/protobuf/json/GFA) index: index features of the graph in a disk-backed key/value store find: use an index to find nodes, edges, kmers, or positions paths: traverse paths in the graph align: local alignment map: global alignment (kmer-driven) stats: metrics describing graph properties join: combine graphs (parallel) concat: combine graphs (serial) ids: id manipulation kmers: generate kmers from a graph sim: simulate reads by walking paths in the graph mod: various transformations of the graph surject: force graph alignments into a linear reference space msga: construct a graph from an assembly of multiple sequences validate: determine if graph is valid filter: filter reads out of an alignment pileup: pileup reads onto graph positions and edges call: call graph positions from a pileup ``` ## グラフの構築 ``` vg construct -r small/x.fa -v small/x.vcf.gz >x.vg ``` ## 可視化、変換 ### vgからgfaへの変換 ``` vg view x.vg >x.gfa ``` ### vgからdotへの変換 ``` vg view -d x.vg >x.dot ``` ### gam形式からjsonへの変換 ``` vg view -a x.gam >x.json ``` ## アライメント ``` vg align -s CTACTGACAGCAGAAGTTTGCTGTGAAGATTAAATTAGGTGATGCTTG x.vg ``` ## パイプを利用してfastaファイル + vcfファイルを入力としてアライメント ``` vg construct -r small/x.fa -v small/x.vcf.gz | vg align -s CTACTGACAGCAGAAGTTTGCTGTGAAGATTAAATTAGGTGATGCTTG - ``` ## マッピング(書きかけ) ### グラフの構築 ``` vg construct -r small/x.fa -v small/x.vcf.gz >x.vg ``` ### xg / gcsa インデックスペアの構築 ``` vg index -x x.xg -g x.gcsa -k 11 x.vg ``` ### ※注意 このインデックスを構築するプロセスにおいて1024 bpを超えるノードが存在するグラフはGCSA2ではインデックスを構築できない。そのままindexを構築しようとすると ``` error:[for_each_gcsa_kmer_position_parallel] Graph contains nodes longer than 1024 bp, which can't be indexed by GCSA2. Preprocess the graph with "vg mod -X 1024 old.vg > new.vg" to divide these nodes. note: node 2 is 10000 bp ``` というエラーメッセージが出るので ``` vg mod -X 1024 old.vg > new.vg ``` を実行してからインデックスの構築を行う。 ### rocksdb backed indexを用いたインデックスの構築(xg / gcsa インデックスの代替) ``` vg index -s -k 11 -d x.vg.index x.vg ``` ### readをindexされたグラフにアライメント ``` vg map -s CTACTGACAGCAGAAGTTTGCTGTGAAGATTAAATTAGGTGATGCTTG -x x.xg -g x.gcsa -k 22 >read.gam ``` ### ## variant検出 ## vgの使用例 ### paired-end readのfastqファイルのアライメント #### gfaをvgに変換 ``` vg view -v -F x.gfa > x.vg ``` #### indexを張る ``` vg index -x x.xg -g x.gcsa -k 11 x.vg ``` #### paired-endのfastqファイルをアライメント ``` vg map -x index.xg -g index.gcsa -f reads1.fq -f reads2.fq > reads.gam ``` # ファイルフォーマットについて ## GFAフォーマット ## gamフォーマット vgでアライメント結果を出力する時に用いられているファイルフォーマット。バイナリファイルでありJSONをプロトコルバッファーでシリアライズしている。例えば ``` $ vg view -a mapped.gam|jq ``` のようにパイプで繋いでJSONとして処理すると ``` { "sequence": "TATGATGGCGGAGGATTCACCGGAATGCTGCTCGCCTACA", "path": { "mapping": [ { "position": { "node_id": 6, "offset": 16, "is_reverse": true }, "edit": [ { "from_length": 16, "to_length": 16 } ], "rank": 1 }, { "position": { "node_id": 5, "is_reverse": true }, "edit": [ { "from_length": 24, "to_length": 24 } ], "rank": 2 } ] }, "name": "e621e9227b40db33_2", "fragment_prev": { "name": "e621e9227b40db33_1" }, "identity": 1 } ``` のような結果が得られる。 ``` sequence:マップされたリード path:マップの状況 mapping position node_id:マップされたノードのID is_reverse:ストランドがマイナスのものが当たっていれば、これが表示される edit from_length: to_length: rank:よくわからない name:リードのID fragment_prev/next:paired-endの相方の情報 name:相方のリードのID identity:リードの何割がリファレンスグラフと一致していたか ``` もしくは ``` vg view -a aln.gam >aln.json ``` としてもアライメント結果をJSONとして見ることが出来る。 元のリニアなゲノム上でのポジションを計算するためには"mapping"内の"position"と"edit"の情報を参考にすれば良いようだ。"node_id"からグラフのどのノードかが分かり、"edit"からノード内でそのポジションが何塩基目かが分かる。これらの情報から注目したいポジションが元のリニアなゲノム内で何塩基だったかを算出することが出来る。"node_id"に関してはvg modを使用した場合は元のGFAと対応していないのでもう一度GFAに変換後確認する必要がある。 #### ※大和田さんのQiitaの記事を参考にしました。 # ドキュメント編集の方針 vg使用者どうしで共同編集できるようにしたい。 vgに関して新しい情報が分かりしだい、説明や新しい機能などを補足していく。 # 参考 ・[vg github](https://github.com/vgteam/vg) ・[vg wiki](https://github.com/vgteam/vg/wiki) ・[vgのスタックバイナリを作る](https://kalab.docbase.io/posts/305038)