植物泛基因组分析之可视化vg构建的Pan-SV初次尝试

vg的github链接

已经编译好的程序，下载就能用

sequenceTubeMap 的github主页

这个是vg团队开发的对vg构建的图泛基因组进行可视化的一个工具，按照页面的帮助文档 A JavaScript module for the visualization of genomic sequence graphs

我个人的理解是JavaScrip就是网页工具

vg团队做了一个现成的

这里也可以选择上传自己构建好的图形泛基因组进行可视化

构建这个Pan-SV需要准备的数据是一个参考基因组和基于参考基因组做变异检测的得到的vcf文件

vg的github主页提供了一个示例数据

参考基因组 tiny.fa

>x
CAAATAAGGCTTGGAAATTTTCTGGAGTTCTATTATATTCCAACTCTCTG

vcf文件

##fileformat=VCFv4.1
##fileDate=20141110
##source=mutatrix population genome simulator
##seed=1415643582
##reference=x.fa
##phasing=true
##commandline=mutatrix --dry-run -s 0.05 -i 0.01 -p 2 x.fa
##INFO=<ID=AC,Number=A,Type=Integer,Description="Alternate allele count">
##INFO=<ID=TYPE,Number=A,Type=String,Description="Type of each allele (snp, ins, del, mnp, complex)">
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples at the site">
##INFO=<ID=NA,Number=1,Type=Integer,Description="Number of alternate alleles">
##INFO=<ID=LEN,Number=A,Type=Integer,Description="Length of each alternate allele">
##INFO=<ID=MICROSAT,Number=0,Type=Flag,Description="Generated at a sequence repeat loci">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 1
x 9 . G A 99 . AC=1;LEN=1;NA=1;NS=1;TYPE=snp GT 1|0
x 10 . C T 99 . AC=2;LEN=1;NA=1;NS=1;TYPE=snp GT 1|1
x 14 . G A 99 . AC=1;LEN=1;NA=1;NS=1;TYPE=snp GT 1|0
x 34 . T A 99 . AC=2;LEN=1;NA=1;NS=1;TYPE=snp GT 1|1
x 39 . T A 99 . AC=1;LEN=1;NA=1;NS=1;TYPE=snp GT 1|0

这个vcf文件里都是SNP，尝试手动添加一个ins和del

##fileformat=VCFv4.1
##fileDate=20141110
##source=mutatrix population genome simulator
##seed=1415643582
##reference=x.fa
##phasing=true
##commandline=mutatrix --dry-run -s 0.05 -i 0.01 -p 2 x.fa
##INFO=<ID=AC,Number=A,Type=Integer,Description="Alternate allele count">
##INFO=<ID=TYPE,Number=A,Type=String,Description="Type of each allele (snp, ins, del, mnp, complex)">
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples at the site">
##INFO=<ID=NA,Number=1,Type=Integer,Description="Number of alternate alleles">
##INFO=<ID=LEN,Number=A,Type=Integer,Description="Length of each alternate allele">
##INFO=<ID=MICROSAT,Number=0,Type=Flag,Description="Generated at a sequence repeat loci">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  1
x       9       .       G       A       99      .       AC=1;LEN=1;NA=1;NS=1;TYPE=snp   GT      1|0
x       10      .       C       T       99      .       AC=2;LEN=1;NA=1;NS=1;TYPE=snp   GT      1|1
x       14      .       G       A       99      .       AC=1;LEN=1;NA=1;NS=1;TYPE=snp   GT      1|0
x       20      .       T       AAGGC   99      .       AC=1;LEN=5;NA=1;NS=1;TYPE=indel GT      1|1
x       28      .       TTC     T       99      .       AC=1;LEN=3;NA=1;NS=1;TYPE=indel GT      1|1
x       34      .       T       A       99      .       AC=2;LEN=1;NA=1;NS=1;TYPE=snp   GT      1|1
x       39      .       T       A       99      .       AC=1;LEN=1;NA=1;NS=1;TYPE=snp   GT      1|0

对vcf文件进行压缩，构建索引

bgzip tiny.vcf
tabix tiny.vcf.gz

构建图形泛基因组

~/biotools/vg/vg construct -v tiny.vcf.gz -r tiny.fa > x.vg

vg文件转为gfa文件

~/biotools/vg/vg view x.vg > x.gfa

用bandage看下这个结果

vg文件转换成xg

~/biotools/vg/vg view -d x.vg > x.xg

这个xg和vg有啥去呗以及为啥要转化暂时还搞不明白

https://vgteam.github.io/sequenceTubeMap/ 这个链接里上传自己的数据写的是要上传xg格式，但是我上传的时候遇到了一个报错

暂时搞不清楚是啥原因了

sequenceTubeMap 也可以在自己的服务器配置，按照github主页的教程

git clone https://github.com/vgteam/sequenceTubeMap.git

把仓库克隆下来，然后进入 sequenceTubeMap

然后是

yarn install

这一步遇到了报错，但是忘记截图了，我按照这个链接

sudo apt remove cmdtest
sudo apt remove yarn
curl -sS https://dl.yarnpkg.com/debian/pubkey.gpg | sudo apt-key add -
echo "deb https://dl.yarnpkg.com/debian/ stable main" | sudo tee /etc/apt/sources.list.d/yarn.list
sudo apt-get update
sudo apt-get install yarn -y

这个需要root权限

然后是

yarn build

这一步又是报错，忘记截图了，大体是Node的版本不对

解决办法的链接也找不到了，最后在普通账户下也没有操作成功，在root账号下成功了

vg需要添加到环境变量，在root账号下是把vg放到了/usr/bin/目录下

启动这个网页

yarn serve

这里默认的端口是3000

在自己本地电脑

ssh -N -L 3000:localhost:3000 root@ip

然后再自己本地浏览器输入

localhost:3000

启动的页面

上传自己刚刚生成的xg文件

如果要展示不同的节点还需要上传gbwt文件

这里怎么从xg文件得到gbwt文件暂时搞不清楚了，可以利用gbz文件得到

命令是

~/biotools/vg/vg autoindex --workflow giraffe -r tiny.fa -v tiny.vcf.gz -p abc
~/biotools/vg/vg gbwt -o abc.gbwt -Z abc.giraffe.gbz

加上gbwt文件以后图就变成了如下这样

但是我加 abc.gbw好像没有起作用，还有好多不明白的地方

参考链接

vg index那一步是可以产生gbwt文件的

~/biotools/vg/vg index -x x.xg -g x.gcsa -G x.gbwt -T x.vg -t 8

这里生成x.xg和x.gbwt

不过这个加载到网站上达不到上面的效果

vg的可视化命令

~/biotools/vg/vg viz -x t.xg --out graph.svg