SNV_counter.sh

SNV_counter.sh

#!/bin/bash

export LC_ALL=C

printf "*** SNVcounter\n
	*** 'samtools, bcftools, and bedtools' must be installed to run this script\n
	*** Prepare 'reference fasta' and 'SNV list' in BED format\n
	*** merge consecutive SNVs using\n
	*** 'sort -k1,1 -k2,2n SNV.bed | mergeBed -c 4 -o collapse'\n
	*** usage : ./SNVcounter Input.Bam Reference.Fasta mergedSNV.bed\n\n"

printf "()()() Processing '$1' \n\n"

dir=$(pwd)
if [ ! -d "$dir" ]; then
mkdir tmp
fi

cd tmp

samtools mpileup -d 10000000 -Buf ../$2 ../$1 | bcftools view -bvcg - > $1.raw.bcf

# to get information on all sites, remove "v" option in bcftools
#samtools mpileup -d 10000000 -Buf ../$2 ../$1 | bcftools view -bcg - > $1.raw.bcf


# "-B" option to report the base composition at every mapped position

bcftools view $1.raw.bcf | grep -v "#" | cut -f 1,2 > $1.id
bcftools view $1.raw.bcf | grep -v "#" | cut -f 1,8 |sed 's/;/\t/g' | grep -Po 'DP=\d+' > $1.dp
bcftools view $1.raw.bcf | grep -v "#" | cut -f 1,8 |sed 's/;/\t/g' | grep -Po 'DP4=\d+,\d+,\d+,\d+' > $1.dp4

paste $1.id $1.dp $1.dp4 | sed 's/ \|=\|,/\t/g' > $1.paste 

awk 'OFS="\t"{if($6+$7+$8+$9 >= 100) print $1,$2-1,$2-1,FILENAME,$4,$6+$7,$8+$9}' $1.paste |sed 's/.trimmed.SE.fq.gz.sam.sorted.bam.paste//g' > $1.bed

intersectBed -wa -wb -loj -a ../$3 -b $1.bed | cut -f 1,2,3,8,9,10,11 | sort -k1,1 -k2,2n | mergeBed -c 4,5,6,7 -o first,mean,mean,mean > $1.intersect

printf "\n()()() Done '$1' !\n\n"

cd ..