推文的代码主要来源于论文
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文件就可以将二代测序数据比对上去对结构变异进行基因分型了
以上只是最简单的流程，没一个步骤都有很多个参数，大部分参数都没有搞明白是什么意思，还需要花时间再看
推文记录的是自己的学习笔记，很可能存在错误，请大家批判着看
