Pipeline

## Each part of this pipeline was run separately
## This will probably not run as-is, I pulled each part of the pipeline into a single script for convenience
#!/bin/bash

base_dir=/scratch/asherrar/vg_test
ref=/home/asherrar/t2t_sequence/v2.0/chm13v2.0.fa
vardir=/scratch/asherrar/sgdp
tmpdir=/scratch/asherrar/vg_tmp

cd $base_dir

variants=""
for i in {1..22} X Y; do variants+="-v $vardir/SGDP.CHM13v2.0.chr$i.recalibrated.snp_indel.pass.vcf.gz "; done

## < Step 1 - Graph Construction > ##

# construct base graph and index
echo "constructing base graph"
vg autoindex -r $ref -p $base_dir/t2t $variants -T $tmpdir -M 256G -V 2
vg convert -p $base_dir/t2t.vg > $base_dir/t2t.xg
vg convert -p $base_dir/t2t.xg > $base_dir/t2t.pg

## < Step 2 - Read Mapping > ##

echo "mapping reads"
vg map -m "long" -d $base_dir/t2t -f /scratch/asherrar/aicardi_longread/PAE47322.fasta > $base_dir/PAE47322.gam

## < Step 3 - Graph Chunking > ##

mkdir -p $base_dir/chunks
vg chunk -x $base_dir/t2t.xg -M -O pg -b $base_dir/chunks/t2t

## < Step 4 - Calling against Chunks > ##

# run for each chunk:

base_dir=/scratch/asherrar/vg_test
mkdir -p $base_dir/chunk_steps
mkdir -p $base_dir/chunk_vcf
chunkdir=$base_dir/chunk_steps

chunkname=chr${current_chunk}
chunk=$base_dir/chunks/t2t_$chunkname.pg

echo "calling variants on $chunkname"

echo " - augmenting graph"
vg augment -v -p -s -t 24 -m 4 -q 5 -Q 5 -A $chunkdir/${chunkname}_PAE47322_augmented.gam $chunk $base_dir/PAE47322.gam > $chunkdir/${chunkname}_PAE47322_augmented.pg
echo " - calculating snarls"
vg snarls $chunkdir/${chunkname}_PAE47322_augmented.pg > $chunkdir/${chunkname}_PAE47322_augmented.snarls
echo " - packing graph"
vg pack -x $chunkdir/${chunkname}_PAE47322_augmented.pg -g $chunkdir/${chunkname}_PAE47322_augmented.gam -o $chunkdir/${chunkname}_PAE47322_augmented.pack
echo " - calling variants"
vg call $chunkdir/${chunkname}_PAE47322_augmented.pg -r $chunkdir/${chunkname}_PAE47322_augmented.snarls -k $chunkdir/${chunkname}_PAE47322_augmented.pack -s PAE47322_reads > $base_dir/chunk_vcf/${chunkname}_PAE47322_readcalls.vcf