bampeTobdg

#! /bin/bash

# Vivek Rai <vivek.rai@roche.com>
# May 7, 2023
#
# Convert paired-end BAM files to BedGraph file.
# Write to the stdout if output file is not specified.

usage() {
	echo "Usage: $0 -i <input.bam> -o <output.bedgraph> -g <genome> -h <help>"
	exit 1
}

error() {
	echo "$@" 1>&2
}

while getopts "i:o:g:h" opt; do
	case $opt in
	i) input=$OPTARG ;;
	o) output=$OPTARG ;;
	g) genome=$OPTARG ;;
	h) usage ;;
	*) usage ;;
	esac
done

if [ -z "$input" ]; then
	usage
fi

if [ ! -f "$input" ]; then
	echo "Input file $input does not exist"
	exit 1
fi

if [ -f "$output" ]; then
	echo "Output file $output already exists"
	exit 1
fi

if [ ! -f "$genome" ]; then
	echo "Genome file $genome does not exist"
	exit 1
fi

# Get the sample name from input file bath by stripping .bam extension
sample=$(basename "$input" .bam)
original_input="$input"

# Check if BAM file is name sorted and offer to sort if not
if [ "$(samtools view -H "$input" | grep -c 'SO:queryname')" -eq 0 ]; then
	error "BAM file is not name sorted. Sorting now..."
	samtools sort -n "$input" -o "$sample".sorted.bam
	input="$sample".sorted.bam
fi

# Convert BAM to BED
error "1/4 Converting BAM to BED..."
bedtools bamtobed -bedpe -i $input >"$sample".bed 2>"$sample.bampeTobdg.log"

# [...]then selecting the 5' and 3' coordinates of the read pair to generate a
# new BED3 file. Filter improper size fragments

error "2/4 Creating fragments BED file..."
awk '$1==$4 && ($6-$2 > 0) {print $0}' "$sample.bed" | cut -f1,2,6 |
	sort -k1,1 -k2,2n -k3,3n >"$sample.fragments.bed"

# # Remove chromosomes that are not in the genome file
# error "3/4 Removing chromosomes that are not in the genome file..."
# awk 'NR==FNR{a[$1];next} ($1 in a)' "$genome" $sample.fragments.bed >$sample.clean.bed

# Convert BED to BedGraph and write to output if available in an idiomatic manner
error "3/4 Converting BED to BedGraph..."
if [ -z "$output" ]; then
	bedtools genomecov -bg -i $sample.fragments.bed -g "$genome"
else
	bedtools genomecov -bg -i $sample.fragments.bed -g "$genome" >"$output"
fi

error "4/4 Cleanup."
rm "$sample".fragments.bed
if [ "$input" != "$original_input" ]; then
    rm "$input"
fi