#!/usr/bin/perl

# Use HOMER for ERa peak calls.
# Make "original" (as is from HOMER) and re-defined 200mer versions. -- UPDATE: Use 200mer only. --


$tool = "/ddn/gs1/home/grimmsa/tools/ucsc/bedToBigBed";
$cs = "/ddn/gs1/shared/fargod/reference_genomes/mm10/mm10ordered.chromSizes";

$chrlist = `cut -f1 $cs`;

$peakdir = "/ddn/gs1/home/grimmsa/sylvia/histonemod_Aug2016/HOMER_peaks/peak_calls";
#@origfiles = ("21d_E2_ERa.homer-peakcalls.factor_F8_e-5.bed", "21d_Veh_ERa.homer-peakcalls.factor_F8_e-5.bed", 
#"adult_E2_ERa.homer-peakcalls.factor_F8_e-5.adultInput.bed", "adult_E2_ERa.homer-peakcalls.factor_F8_e-5.origInput.bed", 
#"adult_Veh_ERa.homer-peakcalls.factor_F8_e-5.adultInput.bed", "adult_Veh_ERa.homer-peakcalls.factor_F8_e-5.origInput.bed");
@origfiles = ("21d_E2_ERa.homer-peakcalls.factor_F12_e-5.bed", "21d_Veh_ERa.homer-peakcalls.factor_F12_e-5.bed",
"adult_E2_ERa.homer-peakcalls.factor_F12_e-5.origInput.bed", "adult_Veh_ERa.homer-peakcalls.factor_F12_e-5.origInput.bed");


foreach $origfile (@origfiles) {
  $bedfile = "$peakdir/$origfile";
#  $bbfile = "./$origfile"; $bbfile =~ s/\.bed/\.bigBed/;
#  unless (-e $bbfile) { print "$tool $bedfile $cs $bbfile &\n"; }

  $newfile = "$origfile"; $newfile =~ s/\.bed/\.200nt\.bed/;
  if ($newfile =~ /origInput/) { $newfile =~ s/\.origInput//; }
  open(IN, "$bedfile");
  open(OUT, ">$newfile.tmp");
  while (<IN>) {
    chomp $_; ($chr, $p0, $p2, @rest) = split/\t/, $_;
    $p1 = $p0+1;
    $midpt = int(($p1+$p2)/2);
    $newp0 = $midpt - 100;
    $newp2 = $midpt + 100;
    print OUT "$chr\t$newp0\t$newp2\n";
  }
  close(IN); close(OUT);
  if (-e $newfile) { system "rm $newfile"; }
  foreach $chr (split/\n/, $chrlist) { system "grep -w $chr $newfile.tmp | sort -k2,2n -k3,3n >> $newfile"; }
  system "rm $newfile.tmp";
  $newbb = $newfile; $newbb =~ s/\.bed/\.bigBed/;
  unless (-e $newbb) { print "$tool $newfile $cs $newbb &\n"; }
}