#!/bin/bash

#Requires bam files and bed file containing coordinates of peaks of interest
#Determines regions within peaks where included fragement coverage best correlates
#with total read count under peak.  Exports regions as bed file.

#required arguements:
#($1): series of bam files
#($2): file names of bam files
#($3): window bed file
#($4): percent of window fragments must cover, default is 0.9
#($5): output coverage file
#($6): output spanning fragment coverage file

#Usage: ATAC_qPCR_regions.bash bam_series name_list qPCR_windows.bed 0.9 combined.o.bed.coverage combined.f9.bed.coverage

bam_series=$1
name_list=$2
qPCR_windows=$3
percent_overlap=${4:-0.9}
combined_o=$5
combined_f9=$6

COUNTER=0
echo $bam_series |tr , '\n' > bam_series.txt
echo $name_list | tr , '\n' > name_list.txt
paste name_list.txt bam_series.txt > bam_name_list.txt
sort -k1 bam_name_list.txt > sorted_bam_name_list.txt
cut -f 2 sorted_bam_name_list.txt > sorted_bam_list.txt

cat sorted_bam_list.txt | while read file
do

#create bed file of fragments for each bam file

samtools sort -n -o - $file | samtools fixmate - - | bamToBed -i - -bedpe | cut -f 1,2,6 > "${file}".includeinsert.bed

#determine window coverage for each bam file
#combined_o: Total reads in peak/window
#combined_f9: Spanning reads in window

if [ $COUNTER -eq 0 ]
then
coverageBed -a $qPCR_windows -b $file | cut -f 1,2,3,4,5 > $combined_o
coverageBed -a $qPCR_windows -b "${file}".includeinsert.bed -f $percent_overlap | cut -f 1,2,3,4,5 > $combined_f9

else
coverageBed -a $qPCR_windows -b $file | cut -f 5 > temp.o.bed.coverage
coverageBed -a $qPCR_windows -b "${file}".includeinsert.bed -f $percent_overlap | cut -f 5 > temp.f9.bed.coverage
paste $combined_o temp.o.bed.coverage > combined.temp.o.bed.coverage
mv combined.temp.o.bed.coverage $combined_o
paste $combined_f9 temp.f9.bed.coverage > combined.temp.f9.bed.coverage
mv combined.temp.f9.bed.coverage $combined_f9

rm temp.o.bed.coverage
rm temp.f9.bed.coverage

fi

rm "${file}.includeinsert.bed"

let COUNTER=COUNTER+1

done


rm *.txt
