vg.md

graph genomeのツールvgの情報のメモ

例えばbubbleが2つあるようなグラフがあるとする

H       VN:Z:1.1
S       1       GATTAGA
S       2       A
S       3       T
S       4       TTTCGAT
S       5       G
S       6       GTT
S       7       AGGATCG
P       chr1    1+,2+,4+,5+,7+  *
P       chr2    1+,3+,4+,6+,7+  *
W       sample1 0       chr1    0       23      >1>2>4>5>7
W       sample2 0       chr1    0       25      >1>3>4>6>7
L       1       +       3       +       0M
L       1       +       2       +       0M
L       2       +       4       +       0M
L       3       +       4       +       0M
L       4       +       6       +       0M
L       4       +       5       +       0M
L       5       +       7       +       0M
L       6       +       7       +       0M

index

• gfa: グラフ表現(テキスト)
• vg: (編集可能な)グラフ表現(protobufのbinary)
• xg: vgの検索用index file(編集不可能) (vgを簡潔表現にしたもの)
• gbwt: haplotypeをnodeのsequenceとして表現してBWTにしたもの
• gg: 各nodeの塩基を持っている
• gbz=gg+gbwt: この2つが揃って初めてグラフが表現できる

gfa/vg/xg/gbzは相互変換できる

# gfa -> vg
vg convert -g toy.gfa > toy.vg
# vg -> xg
vg index toy.vg -x toy.xg
# xg -> vg
vg view toy.xg
# gfa -> gbz
vg gbwt --gbz-format -g toy.gbz -G toy.gfa
# gbwt+gg -> gbz
vg gbwt -o toy.gbwt -g toy.gg -Z toy.gbz
vg convert -b toy.gbwt toy.gg > toy.vg

など

Reference (例えばCHM13:chr8:100-200でqueryできるようにしたい)とHaplotype(各サンプルの配列情報)は区別されているらしい GFAで言うとP (Path) と W (Walk) が違ったりする GFA v1.1 http://gfa-spec.github.io/GFA-spec/GFA1.html

https://github.com/vgteam/vg/wiki/VG-GBWT-Subcommand

extract

グラフである領域だけ取り出したい。 部分グラフを取り出したり、対応するhaplotypeの配列を取り出したりしたい

vg find

-Eをつけないと、本当にchr1:0-10上にあるpathしか残らない(他のpathはnodeが飛び飛びになってしまう)

% vg find -p chr1:0-10 -x toy.vg | vg view -
H       VN:Z:1.1
S       1       GATTAGA
S       4       TTTCGAT
S       2       A
P       chr1    1+,2+,4+        *
P       chr2    1+,4+   *
W       sample1 0       chr1    0       15      >1>2>4
W       sample2 0       chr1    0       14      >1>4
L       1       +       2       +       0M
L       2       +       4       +       0M
% vg find -p chr1:0-10 -x toy.vg -E | vg view -
H       VN:Z:1.1
S       1       GATTAGA
S       4       TTTCGAT
S       2       A
S       3       T
P       chr1    1+,2+,4+        *
P       chr2    1+,3+,4+        *
W       sample1 0       chr1    0       15      >1>2>4
W       sample2 0       chr1    0       15      >1>3>4
L       1       +       2       +       0M
L       1       +       3       +       0M
L       2       +       4       +       0M
L       3       +       4       +       0M

vg chunk

あるノードの周囲c個を部分グラフとして切り出せる -c 0の時はvg findと同じようにpath上にあるnodeしか残らない。 gbwtの時、haplotypeのsample名が消えてしまう？

% vg chunk -x toy.xg -T -c 0 -p chr1:0-10 | vg view -
H       VN:Z:1.1
S       1       GATTAGA
S       4       TTTCGAT
S       2       A
P       chr1    1+,2+,4+        *
P       chr2[0] 1+      *
P       chr2[8] 4+      *
W       sample1 0       chr1    0       15      >1>2>4
W       sample2 0       chr1    0       7       >1
W       sample2 0       chr1    8       15      >4
L       1       +       2       +       0M
L       2       +       4       +       0M
% vg chunk -x toy.xg -G toy.gbwt -T -c 0 -p chr1:0-10 | vg view -
H       VN:Z:1.1
S       1       GATTAGA
S       4       TTTCGAT
S       2       A
P       chr1    1+,2+,4+        *
P       chr2[0] 1+      *
P       chr2[8] 4+      *
P       thread_0        1+,2+   *
P       thread_1        1+      *
W       sample1 0       chr1    0       15      >1>2>4
W       sample2 0       chr1    0       7       >1
W       sample2 0       chr1    8       15      >4
L       1       +       2       +       0M
L       2       +       4       +       0M

vgteam/vg#3258 にも書いてあるように、odgiにしてから切るのが一番楽そうだ

参考

READMEとwiki https://github.com/vgteam/vg https://github.com/vgteam/vg/wiki

vgのアドベントカレンダー https://github.com/genomegraph/vg-advent-calendar-2019