#this script report the commands used to call variants from the 16 Illumina WGS barn swallow individuals.
#for steps 1-3, this github solution was followed: https://github.com/vgteam/vg/issues/3411
#for steps 4-9, the section "Calling the graph by first splitting into components" of this tutorial was followed: https://github.com/vgteam/vg/wiki/Whole-genome-calling-and-genotyping
mkdir tmp
#1. Chop nodes in the graph so they are not more than 256 bp and index the graph
vg mod -t 32 -X 256 graph.vg > graph_mod_chopped.vg
vg index -t 32 -x graph_mod_chopped.xg graph_mod_chopped.vg
#2. Prune the graph with kmer size 45 and index the graph
vg prune -t 32 -k 45 graph_mod_chopped.vg > graph_mod_chopped_pruned.vg
vg index -b ./tmp -t 32 -p -L -x graph_mod_chopped_pruned.xg -g graph_mod_chopped_pruned.gcsa graph_mod_chopped_pruned.vg
vg convert graph.xg -p > graph.pg ## convert the graph (can be vg or xg)
#3. Align the individuals Illumina WGS data
for sample in C24 Cvi Kyo Sha;do
# 比对
vg map -m short \
-t 32 \
-x graph_mod_chopped_pruned.xg \
-g graph_mod_chopped_pruned.gcsa \
-f $datadir/fastq/$sample.illumina.reads_1.fq.gz \
-f $datadir/fastq/$sample.illumina.reads_2.fq.gz > $sample.gam
vg stats -a $sample.gam > $sample.gam.stats
# Augment the subgraph 扩展 graph
vg augment -t 32 -Q 5 -q 5 -m 4 -s -A $sample.aug.gam graph.vg $sample.gam > $sample.aug.vg
# 基因组大，建议分析snarls，节约时间
vg snarls -t 32 $sample.aug.vg > $sample.aug.snarls
# 构建 xg index
vg index -t 32 -L -b tmp -x $sample.aug.xg $sample.aug.vg
# filter secondary and ambiguous read mappings out of the gam
vg filter -t 32 $sample.aug.gam -r 0.90 -fu -m 1 -q 15 -D 999 -x $sample.aug.xg > $sample.filtered.gam
vg gamsort -t 32 -p $sample.filtered.gam -i $sample.filtered.sorted.gam.gai > $sample.filtered.sorted.gam
# Compute the read support, read 覆盖情况
vg pack -t 32 -x $sample.aug.xg -g $sample.filtered.sorted.gam -Q 5 -s 5 -o $sample.aug.pack
# Call variants 变异检测
# vg call -t 32 -a -i -s $sample -k $sample.aug.pack -r $sample.aug.snarls $sample.aug.xg > $sample.aug.vcf
vg call -t 32 -a -k $sample.aug.pack -r $sample.aug.snarls $sample.aug.xg > $sample.aug.vcf
done
