

# Make universe of Peak regions
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/ES_D10/H3K4me1/overlap-peaks.sh-1.0.0/experiment_pooled.sorted.narrowPeak.replicated.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/ES_D10.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/ES_D2/H3K4me1/overlap-peaks.sh-1.0.0/experiment_pooled.sorted.narrowPeak.replicated.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/ES_D2.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/ES_D0/H3K4me1/overlap-peaks.sh-1.0.0/experiment_pooled.sorted.narrowPeak.replicated.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/ES_D0.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/ES_D7/H3K4me1/overlap-peaks.sh-1.0.0/experiment_pooled.sorted.narrowPeak.replicated.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/ES_D7.narrowPeak
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/ES_D5/H3K27ac/overlap-peaks.sh-1.0.0/experiment_pooled.sorted.narrowPeak.replicated.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/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_h3k4me1/ES_D0.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/ES_D2.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/ES_D5.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/ES_D7.narrowPeak /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/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_h3k4me1/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_h3k4me1/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_h3k4me1/overlap_peaks.narrowPeak​ >> /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/unique_peaks.narrowPeak; done
sort-bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/unique_peaks.narrowPeak | cut -f1,2,3 > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/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_h3k4me1/overlap_peaks.narrowPeak​ > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/merge_overlap_peaks.bed
bedops --everything /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/merge_overlap_peaks.bed /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/unique_peaks_sorted.narrowPeak > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/universe_peaks.bed

# Merge peaks that overlapping within 500bp not adjoining
bedtools merge -i /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/universe_peaks.bed -d 500 > /project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/universe_peaks.merge.bed

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

#Identifying genic/intergenic peaks
bedtools intersect -a /Volumes/project/GCRB/Lee_Lab/s163035/Matrix_analysis_PMIT_25842977/ChIP-seq/universe_h3k4me1/universe_peaks.merge.bed -b excluded_3kb_flanking.bed -v > universe_enhancer_H3K4me1.bed

# Get RPKM
./rpkm.py --peaks universe_enhancer_H3K4me1.bed --experiments h3k4me1_list.csv -f H3K4me1
./rpkm.py --peaks universe_enhancer_H3K4me1.bed --experiments h3k4me1_list.csv -f H3K4me1_all --minimum 0

# Graph cutoffs
./cutoff_analysis.py --rpkm H3K4me1_all_filtered_peaks.tsv --color '#008D14' --factor H3K4me1 -l 1
