# Make universe of Peak regions
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/H3K4me3/ES_D10/call-peaks-macs2.sh-2.1.0/ERR207999.fastq.gz.sorted.bam.filtered.no_dups_peaks.narrowPeak> /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D10.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/H3K4me3/ES_D0/call-peaks-macs2.sh-2.1.0/ERR208008.fastq.gz.sorted.bam.filtered.no_dups_peaks.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D0.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/H3K4me3/ES_D2/call-peaks-macs2.sh-2.1.0/ERR208014.fastq.gz.sorted.bam.filtered.no_dups_peaks.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D2.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/H3K4me3/ES_D7/call-peaks-macs2.sh-2.1.0/ERR207987.fastq.gz.sorted.bam.filtered.no_dups_peaks.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D7.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/H3K4me3/ES_D5/call-peaks-macs2.sh-2.1.0/ERR207998.fastq.gz.sorted.bam.filtered.no_dups_peaks.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D5.narrowPeak

# Find overlap of peaks with at least 50% overlap
bedops --everything /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D0.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D2.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D5.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D7.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D10.narrowPeak \
    | bedmap --echo-map --fraction-both 0.5 - \
    | awk '(split($0, a, ";") > 1)' - \
    | sed 's/\;/\n/g' - \
    | sort-bed - \
    | uniq - \
    > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/overlap_peaks.narrowPeak​

# get peaks that are not represented in overlap peaks
for i in $(ls /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/ES_D*narrowPeak); do bedtools intersect -v -wa -f 1.0 -r -a $i -b /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/overlap_peaks.narrowPeak​ >> /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/unique_peaks.narrowPeak; done
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/unique_peaks.narrowPeak | cut -f1,2,3 > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/unique_peaks_sorted.narrowPeak

# Merge Peaks and Union of peaks
bedops -m --range -1 /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/overlap_peaks.narrowPeak​ > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/merge_overlap_peaks.bed
bedops --everything /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/merge_overlap_peaks.bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/unique_peaks_sorted.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/universe_peaks.bed

# Merge peaks that overlapping and adjoining peaks
bedtools merge -i /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/universe_peaks.bed > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/universe_peaks.merge.bed

# Remove extra files
rm /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/overlap_peaks.narrowPeak​
rm /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/unique_peaks.narrowPeak
rm /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/unique_peaks_sorted.narrowPeak
rm /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/merge_overlap_peaks.bed
rm /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/universe_peaks.bed

# Get RPKM
./rpkm.py --peaks /Volumes/project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/universe_peaks.merge.bed --experiments h3k4me3_list.csv -f H3K4me3
./rpkm.py --peaks /Volumes/project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me3/universe_peaks.merge.bed --experiments h3k4me3_list.csv -f H3K4me3_all --minimum 0

# Graph cutoffs
./cutoff_analysis.py --rpkm H3K4me3_all_filtered_peaks.tsv --color '#FF7C21' --factor test -l 1