泛基因组分析流程pggb结合vg流程可以用来做基因组层面的比对并分析snp和indel

推文的代码主要来源于论文

A reference-unbiased super-pangenome of the North American wild grape species (Vitis spp.) reveals genus-wide association with adaptive traits

https://www.biorxiv.org/content/10.1101/2023.06.27.545624v1

葡萄泛基因组.pdf

https://github.com/noecochetel/North_American_Vitis_Pangenome

pggb的github主页

https://github.com/pangenome/pggb

比较详细的文档

https://pggb.readthedocs.io/en/latest/

pggb论文

Building pangenome graphs

https://www.biorxiv.org/content/10.1101/2023.04.05.535718v1

pggb.pdf

接下来用3个拟南芥基因组1号染色体的1:1000 bp试一下基本的命令

首先是两两比对

输入文件压缩且用samtools构建索引

wfmash -t 16 -p 85 -s 10000 -n 1 an1.fasta.gz c24.fasta.gz > an1_c24.paf
wfmash -t 16 -p 85 -s 10000 -n 1 an1.fasta.gz ler.fasta.gz > an1_ler.paf
wfmash -t 16 -p 85 -s 10000 -n 1 c24.fasta.gz ler.fasta.gz > c24_ler.paf

cat *.paf > all.on.all.paf

paf文件结合fasta文件构建gfa文件

zcat an1_chr1.fa.gz c24_chr1.fa.gz ler_chr1.fa.gz > all.on.all.fasta
time seqwish -p all.on.all.paf -s all.on.all.fasta -g all.on.all.chr1.gfa -t 16 -P

接下来是用 smoothxg 命令对获得的gfa文件进行处理，这一步命令的参数很多，暂时没看明白是什么意思

先跳过这一个步骤

vg stats 可以统计gfa文件的基本信息

vg stats -z all.on.all.chr1.gfa

输出 nodes和 edges

gfa文件可以转化成vg文件

time vg view --gfa-in all.on.all.chr1.gfa --vg --threads 16 > all.on.all.chr1.vg

vg文件可以index成xg文件

vg index all.on.all.chr1.vg -x all.on.all.chr1.xg --threads 36

xg文件可以转换成vcf文件

vg deconstruct -p an1_chr1 -t 8 -v -a -e all.on.all.chr1.xg > all.on.all.chr1.vcf

-p 参数就是参考序列的序列id

可以使用bcftools对vcf文件进行基本过滤

比如只保留插入变异

bcftools filter --include 'strlen(REF)<strlen(ALT)' --output all.on.all.ins.vcf --output-type v --threads 4 all.on.all.vcf

这个vcf文件还可以结合参考的fasta 重新构建vg文件

vg construct -r an1_chr1.fa -v all.on.all.ins.vcf.gz -t 8 -m 32 > ins.vg

有了这个vg文件就可以将二代测序数据比对上去对结构变异进行基因分型了

以上只是最简单的流程，没一个步骤都有很多个参数，大部分参数都没有搞明白是什么意思，还需要花时间再看

image.png

推文记录的是自己的学习笔记，很可能存在错误，请大家批判着看