vcfcmp.sc

#!/bin/bash

usage () {
	echo "Usage: $0 ref.vcf var1.vcf var2.vcf ..."
	echo " e.g.: $0 WT.vcf MU.vcf"
	exit 1
}

if [ $# -lt 2 ] ; then
	usage
fi

if [ $1 == "-h" ] ; then
	usage
fi

ref=$1
shift

if [ ! -f $ref ]; then
	echo "Reference VCF file ($ref) does not exist."
	usage
fi

# ref
if [ ! -f $ref.gz ]; then
	bgzip $ref
fi
if [ ! -f $ref.gz.csi ]; then
	bcftools index $ref.gz
fi

mkdir -p ${PWD}/sector

# var
for var in $@; do

	if [ ! -f $var ]; then
		echo "Mutant VCF file ($var) does not exist."
		usage
	fi

	# Use bcftools to subtract background
	echo "Removing background variants in $var ..."
	if [ ! -f $var.gz ]; then
		bgzip $var
	fi
	if [ ! -f $var.gz.csi ]; then
		bcftools index $var.gz
	fi
	bcftools isec $ref.gz $var.gz -p filtered

	# Restore original VCF file
	gunzip $var.gz
	rm -rf $var.gz.csi

	# Use filter_vcf.pl to filter and annotate VCF file
	filter_vcf.sc.pl filtered/0001.vcf

	# Prepare out vcf
	vcf=`basename $var`
	mv ${PWD}/filtered/0001.filtered.ann.vcf ${PWD}/sector/$vcf
	rm -rf filtered

done

# Restore ref.vcf
gunzip $ref.gz
rm -rf $ref.gz.csi

echo "Sector VCF files are available at ${PWD}/sector/"
echo "Done"