#!/usr/bin/perl
use strict;
use warnings;
use File::Path;

# Copyright (c) 2011 Swaraj Basu (swaraj.basu@szn.it)
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ALSO, IT WOULD BE NICE IF YOU LET ME KNOW YOU USED IT.

# TAKE A 12 COLUMN BED FILE AND EXTRACT THE COORDINATES FOR ALL EXONS OF EACH GENE. PRINT THE RESULTS AS 6 COLUMN BED FILE "CHR	START	END GENENAME SCORE STRAND". USE THE EXON FILE TO GET EXON COUNTS AGAINST A SET OF BAM FILES IN A GIVEN DIRECTORY. GIVEN A DIRECTORY OF FILES OF EXON COUNTS ADD TO THE ORIGINAL BED FILE SUM OF ALL EXON COUNTS FOR EACH GENE IN EACH SAMPLE. THE SCRIPT REQUIRES THAT BEDTOOLS AND SAMTOOLS BE PREINSTALLED IN THE SYSTEM. THE SCRIPT TAKES IN PAIRED END READS, TO CHANGE MODIFY LINE 102.

#INPUT FILES
my $bed = "filename.bed";#BED FILE WITH FEATURES FOR WHOM COUNTS ARE REQUIRED.
my $dirBAM = "bam";#DIRECTORY OF BAM FILES (VARIOUS SAMPLES)
my $dirCNT = "tmpcount";
my $exonFILE = $bed.".exon";
my $bedOUT = $bed.".count";


&bed_to_exon($bed, $exonFILE);
&make_exon_count($dirBAM, $exonFILE, $dirCNT);
&sum_exon_count($bed, $dirCNT, $bedOUT);
#rmtree($dirCNT);
#unlink($exonFILE);


#CONVERT BED FILE INTO 12 COLUMN BED FOR EACH EXON OF EACH FEATURE
sub bed_to_exon{
  my $bed = shift;
  my $exonFILE = shift;
  open(FH, "<$bed");
  open(FV, ">>$exonFILE");

  while(<FH>){
    my $line = $_;
    chomp($line);
    my @elements = split("\t", $line);
    my $btype = @elements;
    my $fchr = $elements[0];
    my $fstart = $elements[1];
    my $fend = $elements[2];
    my $fname = $fchr.$fstart.$fend.$elements[3];
    my $fscore = $elements[4];
    my $fstrand = $elements[5];
    my @blsiz = split("\,", $elements[10]);
    my @blstr = split("\,", $elements[11]);
    my $length = @blsiz;
  
    for(my $i = 0; $i < $length; $i++){
      my $bstart = $blstr[$i];
      my $bsize = $blsiz[$i];
      my $estart = $fstart + $bstart;
      my $eend = $fstart + $bstart + $bsize;
      
      #COVERAGE BED WORKS PROPERLY WITH 12 COLUMN BED FILES ONLY
      print FV "$fchr\t$estart\t$eend\t$fname\t$fscore\t$fstrand\t$estart\t$eend\t0,0,0\t1\t$bsize\,\t$bstart\,\n";
    }
  }
  close(FH);
  close(FV);
}

#CHECK THE COUNT OF EACH EXON IN EACH BAM. 
sub make_exon_count{
  my $dirBAM = shift;
  my $exonFILE = shift;
  my $dirCNT = shift;  

  #ACCESS DIRECTORY OF BAM FILES
  opendir(DIR, $dirBAM) || die ( "Cannot open dir!!!" );
  my @bamDIR = readdir(DIR);

  if(-d $dirCNT){
    die "A directory from a previous run exists, please remove or rename it.";
  }
  mkdir($dirCNT) || die "No write permissions to create count directory\n";

  foreach my $bam (@bamDIR){
    next unless $bam =~m/^.*\.bam$/;
    my $bamPATH = $dirBAM."/".$bam;
    my $countOUT = $bam.".zcountz";
    my $countPATH = $dirCNT."/".$countOUT;
    #THE COMMAND CURRENTLY ASSUMES READS TO BE PAIRED END AND FILTERS BY A QUALITY SCORE OF 30. CHANGE THE PARAMETERS. -f is get paried mapped reads. Add -q if quality filtering is needed.
    system("samtools view -u -f 0x2 $bamPATH|coverageBed -split -abam stdin -b $exonFILE|awk \'BEGIN{FS=\"\t\"}\;\{print \$4\"\t\"\$13\}\' \> $countPATH");
  }
}


#GIVEN A DIRECTORY OF FILES OF EXON COUNTS FROM THE SCRIPT  bed_to_exons.pl ALONG WITH THE ORIGINAL BED FILE THIS SCRIPT ADDS TO THE ORIGINAL BED FILE SUM OF ALL EXON COUNTS FOR EACH GENE IN EACH SAMPLE.
sub sum_exon_count {
  my $bed = shift;
  my $dirCNT = shift;
  my $bedOUT = shift;

  open(FH, "<$bed");
  open(FV, ">>$bedOUT");
  my @bed = <FH>;
  close(FH);

  opendir(DIR, $dirCNT) || die ( "Cannot open dir!!!" );
  my @files = readdir(DIR);
  close(DIR);

  my $mahahash;

  foreach my $file (sort @files){
    next if(-B $file);#AVOID BINARY FILES
    my $filePATH = $dirCNT."/".$file;
    
    open(FC, "<$filePATH")|| die ( "Cannot open $file!!!" );
    my @cfile = <FC>;
    close(FC);
    foreach my $line (@cfile){
      chomp($line);
      my @elements = split("\t", $line);
      my $name = $elements[0];
      my $count = $elements[1];
      $mahahash->{$file}->{$name} += $count;
    }
  }

  my $namelist = join ("\t", keys(%{$mahahash}));
  print "$namelist\n";
  print FV "Chromosome\tStart\tEnd\tName\tScore\tStrand\tthickStart\tthickEnd\titemRgb\tblockCount\tblockSizes\tblockStarts\t$namelist\n";
  
  foreach my $line(@bed){
    chomp($line);
    print FV $line."\t";
    my @elements = split ("\t", $line);
    my $fchr = $elements[0];
    my $fstart = $elements[1];
    my $fend = $elements[2];
    my $name = $fchr.$fstart.$fend.$elements[3];
    my $length = keys(%{$mahahash});
    my $i = 1;
    foreach my $key (keys %{$mahahash}){
      
      my $count = $mahahash->{$key}->{$name};
      if($i == $length){
        print FV $count."\n";
      }
      else{
        print FV $count."\t";
      }
      $i++;
    }
  }
close(FH);
close(FV);
}