#!/bin/bash

# Script for quantifying enhancer activity
# 'mm10_chrom_lengths.bed' & 'gencode.vM18.annotation.boundaries.uniq.bed' supplied in repos
# 'gencode.vM18.annotation.boundaries.uniq.bed' can also be generated by running collapse_gene_boundaries.sbatch
# 'picc_BMDM_WT_IFNG_2h_H3K27ac.sorted.bam' generated using scripts in ChIPseq_BMDM repos

# Example usage:
# sbatch find_neighborhoods.sbatch

# General settings
#SBATCH -p short
#SBATCH -N 1
#SBATCH -n 1
#SBATCH --time=6:00:00
#SBATCH --mem=16G

# Job name and output
#SBATCH -J find_neighborhoods
#SBATCH -o /Users/%u/slurmOut/slurm-%j.out
#SBATCH -e /Users/%u/slurmErr/slurm-%j.err

# Load modules
module purge
module load gcc/7.1.0 bedtools/2.28.0 python/3.7.4

# Sort candidateRegions.bed
pwd; hostname; date

echo $(date +"[%b %d %H:%M:%S] Sorting candidateRegions.bed...")

bedtools sort -i plat_BMDM_IFNG_2h_ATAC_R1_peaks.narrowPeak.candidateRegions.bed \
> plat_BMDM_IFNG_2h_ATAC_R1_peaks.narrowPeak.candidateRegions.sorted.bed

# Call candidate regions
echo $(date +"[%b %d %H:%M:%S] Calling candidate regions...")

python ABC-Enhancer-Gene-Prediction-0.2.2/src/run.neighborhoods.py \
--candidate_enhancer_regions plat_BMDM_IFNG_2h_ATAC_R1_peaks.narrowPeak.candidateRegions.sorted.bed \
--genes gencode.vM18.annotation.boundaries.uniq.bed \
--H3K27ac picc_BMDM_WT_IFNG_2h_H3K27ac.sorted.bam \
--ATAC plat_BMDM_IFNG_2h_ATAC_R1.sorted.dedup.bam,plat_BMDM_IFNG_2h_ATAC_R2.sorted.dedup.bam,plat_BMDM_IFNG_2h_ATAC_R3.sorted.dedup.bam \
--chrom_sizes mm10_chrom_lengths.bed \
--outdir ./

echo $(date +"[%b %d %H:%M:%S] Done")