#!/usr/bin/env perl

##
# RF PeakCall
# RNA Framework [http://www.rnaframework.com]
#
# Author:  Danny Incarnato (dincarnato[at]rnaframework.com)
# Summary: Performs peak calling of RNA immunoprecipitation (IP) experiments
#
# This program is free software, and can be redistribute  and/or modified
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# any later version.
#
# Please see <http://www.gnu.org/licenses/> for more informations.
##

use strict;
use Config;
use File::Basename;
use FindBin qw($Bin);
use Getopt::Long qw(:config no_ignore_case);
use threads;
use threads::shared;

use lib $Bin . "/lib";

use Core::Mathematics qw(:all);
use Core::Statistics;
use Core::Utils;
use Data::IO;
use Data::Sequence::Utils;
use Graphics::Chart::Linedot;
use Graphics::Image;
use RF::Data::IO::RC;
use RF::Utils;
use Term::Constants qw(:screen);
use Term::Progress;

use constant SIZE5UTR => 75;
use constant SIZECDS  => 650;
use constant SIZE3UTR => 275;
use constant SIZERNA  => 1000;

$|++;

die "\n  [!] Error: This program requires ithreads." .
    "\n             Please recompile Perl with ithreads support and try again\n\n" unless(defined $Config{useithreads});

my ($help, $output, $overwrite, $win, $tmpDir,
    $pvalue, $pseudocount, $meancov, $mediancov,
    $control, $ip, $index, $enrich, $ctrlIndex,
    $ipIndex, $threads, $iprc, $ctrlrc, $decimals,
    $mergedist, $img, $metaplot, $refine, $summit, 
    $relaxed, $offset, $metaonlyenriched, $whitelist, 
    $metacoding, $minorflen, $altstart, $anystart, 
    $gencode, $codsize, $orfFile, $R, @pool, @imeans, 
    %orfs);

my $ctrlnreads : shared;
my $ipnreads : shared;
my $progressBar : shared;
my @qcaller : shared;
my @qresults : shared;
my @metaip : shared;
my @metactrl : shared;
my @metapeaks : shared;
my @metaipcod : shared;
my @metactrlcod : shared;
my @metapeakscod : shared;
my %results : shared;
%results = ( cov      => 0,
             incov    => 0,
             diffcseq => 0,
             nocid    => 0,
             pass     => 0,
             peaks    => 0,
             codpeaks => 0 );

do {

    local $SIG{__WARN__} = sub { };

    GetOptions( "h|help"                => \$help,
                "o|output=s"            => \$output,
                "ow|overwrite"          => \$overwrite,
                "g|img"                 => \$img,
                "mp|meta-plot"          => \$metaplot,
                "eo|plot-enriched-only" => \$metaonlyenriched,
                "w|window=i"            => \$win,
                "l|whitelist=s"         => \$whitelist,
                "mc|mean-coverage=s"    => \$meancov,
                "ec|median-coverage=s"  => \$mediancov,
                "c|control=s"           => \$control,
                "I|IP=s"                => \$ip,
                "i|index=s"             => \$index,
                "D|decimals=i"          => \$decimals,
                "p|processors=i"        => \$threads,
                "f|offset=i"            => \$offset,
                "e|enrichment=s"        => \$enrich,
                "v|p-value=s"           => \$pvalue,
                "r|refine"              => \$refine,
                "x|relaxed"             => \$relaxed,
                "s|summit"              => \$summit,
                "of|orf-file=s"         => \$orfFile,
                "pc|pseudocount=s"      => \$pseudocount,
                "md|merge-distance=i"   => \$mergedist,
                "mcp|meta-coding-plot"  => \$metacoding,
                "mo|min-orf-len=i"      => \$minorflen,
                "als|alt-start"         => \$altstart,
                "ans|any-start"         => \$anystart,
                "gc|genetic-code"       => \$gencode ) or help(1);

};

help() if ($help);

# Default
$codsize = SIZE5UTR + SIZECDS + SIZE3UTR;
$threads ||= 1;
$win ||= 150;
$offset ||= round($win / 2);
$enrich ||= 3;
$mergedist //= 50;
$pvalue //= 0.05;
$pseudocount ||= 1;
$decimals ||= 3;
$meancov //= 0;
$mediancov //= 0;
$gencode ||= 1;
$minorflen ||= 50;
$R = checkRinstall($R) if ($img || $metacoding || $metaplot);

die "\n  [!] Error: No IP sample RC file provided\n\n" if (!defined $ip);
die "\n  [!] Error: Provided IP sample RC file does not exist\n\n" if (!-e $ip);
die "\n  [!] Error: Provided RCI index does not exist\n\n" if (defined $index && !-e $index);
die "\n  [!] Error: Provided ORF file does not exist\n\n" if (defined $orfFile && !-e $orfFile);
die "\n  [!] Error: Provided control sample RC file does not exist\n\n" if (defined $control && !-e $control);
die "\n  [!] Error: Number of processors must be an integer greater than 0\n\n" if ($threads < 1);
die "\n  [!] Error: Decimals value must be an integer comprised between 1 and 10\n\n" if ($decimals < 1 || $decimals > 10);
die "\n  [!] Error: Window's size must be an integer >= 10\n\n" if ($win < 10);
die "\n  [!] Error: Offset cannot exceed window's size\n\n" if ($offset > $win);
die "\n  [!] Error: P-value must be numeric\n\n" if (!isnumeric($pvalue));
die "\n  [!] Error: P-value must be comprised between 0 and 1\n\n" if (!ispositive($pvalue) || $pvalue > 1);
die "\n  [!] Error: Pseudocount value must be numeric\n\n" if (!isnumeric($pseudocount));
die "\n  [!] Error: Pseudocount value must be > 0\n\n" if (!$pseudocount);
die "\n  [!] Error: Enrichment value must be numeric\n\n" if (!isnumeric($enrich));
die "\n  [!] Error: Enrichment value must be > 0\n\n" if (!$enrich);
die "\n  [!] Error: Mean coverage value must be numeric\n\n" if (!isnumeric($meancov));
die "\n  [!] Error: Mean coverage value must be >= 0\n\n" if (!ispositive($meancov));
die "\n  [!] Error: Median coverage value must be numeric\n\n" if (!isnumeric($mediancov));
die "\n  [!] Error: Median coverage value must be >= 0\n\n" if (!ispositive($mediancov));
die "\n  [!] Error: Parameter -x requires peaks refinement (-r) to be active\n\n" if ($relaxed && !$refine);
die "\n  [!] Error: Parameter -eo requires meta-plot generation (-mp) to be active\n\n" if ($metaonlyenriched && (!$metaplot && !$metacoding));
die "\n  [!] Error: Provided whitelist file does not exist\n\n" if (defined $whitelist && !-e $whitelist);

warn "\n  [!] Warning: No control sample RC file provided. Background will be estimated on IP sample only.\n" if (!defined $control);

if (!defined $output) {

    my ($cid, $iid);
    $cid = fileparse($control, qr/\.[^.]*/) if (defined $control);
    $iid = fileparse($ip, qr/\.[^.]*/);

    $output = (defined $control ? $iid . "_vs_" . $cid : $iid) . "/";

}
else { $output =~ s/\/?$/\//; }

$SIG{__DIE__} = \&cleanup;
$tmpDir = "${output}tmp";

print "\n[+] Making output directory...";

if (-e $output) {

    if ($overwrite) {

        my $error = rmtree($output);

        die "\n\n  [!] Error: " . $error . "\n\n" if ($error);

    }
    else { die "\n\n  [!] Error: Output directory already exists." .
               "\n      Please use -ow (or --overwrite) to overwrite output directory\n\n"; }

}

if (my $error = mktree($tmpDir, "755")) { die "\n  [!] Error: Unable to create output directory ($error)\n\n"; }

mktree($output . "plots", "755") if ($metaplot || $metacoding);
mktree($output . "plots/transcripts", "755") if ($img);

# In case no index has been provided, we rebuild the index before generating the working threads
# The new index will be generated in the same path of the rc file, with extension .rci
$ctrlIndex = defined $index ? $index : (defined $control && -e "$control.rci" ? "$control.rci" : undef);
$ipIndex = defined $index ? $index : (-e "$ip.rci" ? "$ip.rci" : undef);
$ctrlrc = RF::Data::IO::RC->new( file       => $control,
                                 index      => $ctrlIndex,
                                 buildindex => defined $ctrlIndex ? 1 : 0 ) if ($control);
$iprc = RF::Data::IO::RC->new( file       => $ip,
                               index      => $ipIndex,
                               buildindex => defined $ipIndex ? 1 : 0 );

print "\n[+] Loading transcript IDs... ";

@qcaller = $iprc->ids();

if ($whitelist) {

    my %whitelist;

    open(my $fh, "<", $whitelist) or die "\n\n  [!] Error: Unable to read from whitelist file ($!)\n\n";
    while(<$fh>) {

        chomp();
        $whitelist{$_} = 1;

    }
    close($fh);

    @qcaller = grep { exists $whitelist{$_} } @qcaller;

    if (!@qcaller) { die "\n\n  [!] Error: No transcript matched provided whitelist\n\n"; }

}

print scalar(@qcaller) . " transcripts loaded.";

if ($orfFile) {

    print "\n[+] Importing ORF coordinates...";

    open(my $fh, "<", $orfFile) or die "\n\n  [!] Error: Unable to read from ORF file ($!)\n\n";
    while(<$fh>) {

        chomp();

        my @row = split /\t/;

        next if (@row < 3);

        if (my $length = $iprc->length($row[0])) {

            if (isint(@row[1, 2]) && ispositive(@row[1, 2]) && $row[1] < $row[2]) {

                next if ($row[2] > $length);

                $orfs{$row[0]} = [$row[1], $row[2] - 1];

            }
            else { die "\n\n  [!] Error: Malformed ORF file\n\n"; }

        }

    }
    close($fh);

    print " " . scalar(keys %orfs) . " imported.";

}

$progressBar = shared_clone(Term::Progress->new( max     => scalar(@qcaller),
                                                 width   => 50,
                                                 colored => 1 ));

print "\n[+] Calling peaks...\n\n";

$progressBar->init();
@pool = map{ threads->create(\&call) } 1 .. $threads;
$_->join() for(@pool);

print "\n";

if ($results{peaks}) {

    print "\n[+] Sorting peaks...";

    # Sort lexicographically by transcript ID and numerically by start position
    @qresults = sort {$a->[0] cmp $b->[0] ||
                      $a->[1] <=> $b->[1]} @qresults;

    print "\n[+] Writing output BED file...";

    open(my $wh, ">", $output  . "/peaks.bed") or die "\n\n  [!] Error: Unable to write peaks BED file (" . $! . ")\n\n";
    select((select($wh), $|=1)[0]);
    print $wh join("\t", @{$_}[0 .. 4]) . "\n" for (@qresults);
    close($wh);

    if ($summit) {

        open(my $wh, ">", $output  . "/summits.bed") or die "\n\n  [!] Error: Unable to write summits BED file (" . $! . ")\n\n";
        select((select($wh), $|=1)[0]);
        print $wh join("\t", @{$_}[0, 5, 6]) . "\n" for (@qresults);
        close($wh);

    }

}

if ($metaplot &&
    (($metaonlyenriched && @qresults) ||
     !$metaonlyenriched)) {

    my (@icov, @ccov, @pcov);
    @icov = map { ($_ / $results{cov}) / ($metaonlyenriched ? $ipnreads : $iprc->mappedreads()) * (10 ** 6) } @metaip;
    @ccov = map { ($_ / $results{cov}) / ($metaonlyenriched ? $ctrlnreads : $ctrlrc->mappedreads()) * (10 ** 6) } @metactrl if ($control);
    @pcov = map { $_ / @qresults * 100 } @metapeaks if (@qresults);

    plotcoverage($output . "plots/meta-gene.pdf", \@icov, \@ccov, \@pcov);

}

if ($metacoding &&
    (($metaonlyenriched && @qresults) ||
     !$metaonlyenriched)) {

    my (@icov, @ccov, @pcov);
    @icov = map { ($_ / $results{cov}) / (($metaonlyenriched ? $ipnreads : $iprc->mappedreads()) || 1) * (10 ** 6) } @metaipcod;
    @ccov = map { ($_ / $results{cov}) / ($metaonlyenriched ? $ctrlnreads : $ctrlrc->mappedreads()) * (10 ** 6) } @metactrlcod if ($control);
    @pcov = map { $_ / $results{codpeaks} * 100 } @metapeakscod if ($results{codpeaks});

    plotcoverage($output . "plots/meta-coding-gene.pdf", \@icov, \@ccov, \@pcov);

}

cleanup();

print "\n[+] Calling statistics:\n" .
      "\n  [*] Covered transcripts:   " . $results{cov} . " total" .
      "\n                             " . $results{pass} . " transcripts with detected enrichment" .
      "\n                             " . $results{peaks} . " called peaks" .
      "\n  [*] Discarded transcripts: " . ($results{incov} + $results{diffcseq} + $results{nocov} + $results{nocid}) . " total" .
      "\n                             " . $results{incov} . " insufficient coverage" .
      "\n                             " . $results{diffcseq} . " mismatch between IP and control sample sequence" .
      "\n                             " . $results{nocid} . " absent in control sample reference";

print "\n\n[+] All done.\n\n";

sub call {

    my ($crc, $irc, $attributes, $scale,
        @cmeta, @imeta, @cmetacod, @imetacod,
        @pmeta, @pmetacod);

    $crc = RF::Data::IO::RC->new( file  => $control,
                                  index => $ctrlIndex ) if (defined $control);

    $irc = RF::Data::IO::RC->new( file  => $ip,
                                  index => $ipIndex );

    while (1) {

        my ($id, $initial);

        { lock(@qcaller);
          $id = shift(@qcaller) if (@qcaller);

          if ($id && $img) {

            $initial = uc(substr($id, 0, 1));
            mktree($output . "plots/transcripts/" . $initial, "755") if (!-e $output . "plots/transcripts/" . $initial);

          } }

        last unless($id);

        { lock($progressBar);
          $progressBar->update(1); }

        my ($ientry, $sequence, $imedian, $cmedian,
            $last, $length, $ssummit, $esummit,
            $cnreads, $orfStart, $orfEnd, @icov,
            @ccov, @pvalues, @wins, @enrich,
            @peaks);
        $ientry = $irc->read($id);
        $imedian = $ientry->mediancoverage();
        $sequence = $ientry->sequence();
        $length = length($sequence);
        @icov = $ientry->coverage();

        if (defined $control) {

            if (my $centry = $crc->read($id)) {

                $cmedian = $centry->mediancoverage();

                if ($centry->sequence() ne $sequence) {

                    lock(%results);
                    $results{diffcseq}++;

                    next;

                }

                if ($centry->meancoverage() < $meancov ||
                    $cmedian < $mediancov) {

                    lock(%results);
                    $results{incov}++;

                    next;

                }

                $cnreads = $centry->readscount();
                @ccov = $centry->coverage();

            }
            else {

                lock(%results);
                $results{nocid}++;

                next;

            }

        }

        if ($metacoding) {

            if ($orfFile) { if (exists $orfs{$id}) { ($orfStart, $orfEnd) = @{$orfs{$id}}; } }
            else {

                my ($orf);
                ($orf, $orfStart) = longestorf($sequence, { gencode     => $gencode,
                                                            altstart    => $altstart,
                                                            ignorestart => $anystart,
                                                            minlength   => $minorflen });

                $orfEnd = $orfStart + length($orf) - 1 if (defined $orf);

            }

        }

        # Calculates per-base enrichment
        @enrich = $control ? map { logarithm((($icov[$_] + $pseudocount) / ($imedian + $pseudocount)) / (($ccov[$_] + $pseudocount) / ($cmedian + $pseudocount)), 2) } 0 .. $length - 1 :
                             map { logarithm((($icov[$_] + $pseudocount) / ($imedian + $pseudocount)), 2) } 0 .. $length - 1;

        if (!$metaonlyenriched) {

            if ($metaplot) {

                my (@metacoverage);

                if ($control) {

                    @metacoverage = metacoverage(\@ccov, SIZERNA);
                    @cmeta = map { $cmeta[$_] + $metacoverage[$_] } 0 .. SIZERNA - 1;

                    undef(@metacoverage);

                }

                @metacoverage = metacoverage(\@icov, SIZERNA);
                @imeta = map { $imeta[$_] + $metacoverage[$_] } 0 .. SIZERNA - 1;

            }

            if ($metacoding && $orfEnd) {

                my (@cmeta5, @cmeta3, @cmetac, @imeta5,
                    @imeta3, @imetac, @metacod);
                @cmeta5 = (0) x SIZE5UTR;
                @cmeta3 = (0) x SIZE3UTR;
                @imeta5 = @cmeta5;
                @imeta3 = @cmeta3;

                if ($orfStart > 0) {

                    @cmeta5 = metacoverage([@ccov[0 .. $orfStart - 1]], SIZE5UTR) if ($control);
                    @imeta5 = metacoverage([@icov[0 .. $orfStart - 1]], SIZE5UTR);

                }

                @cmetac = metacoverage([@ccov[$orfStart .. $orfEnd]], SIZECDS) if ($control);
                @imetac = metacoverage([@icov[$orfStart .. $orfEnd]], SIZECDS);

                if ($orfEnd < $length - 1) {

                    @cmeta3 = metacoverage([@ccov[$orfEnd + 1 .. $length - 1]], SIZE3UTR) if ($control);
                    @imeta3 = metacoverage([@icov[$orfEnd + 1 .. $length - 1]], SIZE3UTR);

                }

                if ($control) {

                    @metacod = (@cmeta5, @cmetac, @cmeta3);
                    $cmetacod[$_] += ($metacod[$_] ) for (0 .. $#metacod);

                    undef(@metacod);

                }

                @metacod = (@imeta5, @imetac, @imeta3);
                $imetacod[$_] += ($metacod[$_]) for (0 .. $#metacod);

            }

        }

        if ($img) {

            my ($icov, $ccov);
            @{$icov} = map { $_ / $irc->mappedreads() * (10 ** 6) } @icov;
            @{$ccov} = map { $_ / $crc->mappedreads() * (10 ** 6) } @ccov if ($control);

            plotcoverage($output . "/plots/transcripts/" . $initial . "/" . $id . ".pdf", $icov, $ccov);

        }

        for (my $i = 0; $i < $#icov; $i += $offset) {

            my ($start, $end, $score);
            $start = $i;

            if ($start + $win > $#icov) {

                $start = $#icov - $win >= 0 ? $#icov - $win : 0;
                $end = $#icov;

            }
            else { $end = $start + $win; }

            $score = mean(@enrich[$start .. $end]);

            push(@wins, { coords => [$start, $end],
                          score  => [$score],
                          pvalue => undef });

        }

        { lock(%results);
          $results{cov}++; }

        @imeans = map { mean(@icov[$_->{coords}->[0] .. $_->{coords}->[1]]) } @wins;

        foreach my $win (@wins) {

            my $p = defined $control ? fisher(round(mean(@icov[$win->{coords}->[0] .. $win->{coords}->[1]])), round($imedian), round(mean(@ccov[$win->{coords}->[0] .. $win->{coords}->[1]])), round($cmedian)) :
                                       fisher(round(mean(@icov[$win->{coords}->[0] .. $win->{coords}->[1]])), round($imedian), round(mean(@imeans)), round($imedian));
            $win->{pvalue} = $p;

        }

        @pvalues = grep { isnumeric($_) } (map { $_->{pvalue} } @wins);

        next if (!@pvalues);

        @pvalues = padjust(\@pvalues); # Benjamini-Hochberg adjusted p-values

        foreach my $win (@wins) { $win->{pvalue} = [ shift(@pvalues) ] if (isnumeric($win->{pvalue})); }    # Update p-values with BH adjusted p-values

        # Delete peaks whose corrected p-value is greater than cutoff
        # or whose enrichment is lower than cutoff
        for (my $i = 0; $i < @wins; $i++) {

            if (!isnumeric($wins[$i]->{pvalue}->[0]) ||
                $wins[$i]->{pvalue}->[0] >= $pvalue ||
                $wins[$i]->{score}->[0] < $enrich) {

                splice(@wins, $i, 1);

                $i--;

            }

        }

        next if (!@wins);

        if ($refine) {

            $_ = refinepeak( peak       => $_,
                             imedian    => $imedian,
                             cmedian    => $cmedian,
                             length     => $length,
                             icov       => \@icov,
                             ccov       => \@ccov,
                             enrichment => \@enrich ) for (@wins);

        }

        $last = shift(@wins);

        while(my $current = shift(@wins)) {

            if (intersect($current->{coords}, [$last->{coords}->[0], $last->{coords}->[1] + $mergedist])) {

                $last->{coords}->[1] = $current->{coords}->[1];
                push(@{$last->{score}}, $current->{score}->[0]);
                push(@{$last->{pvalue}}, $current->{pvalue}->[0]);

            }
            else {

                push(@peaks, $last);
                $last = $current;

            }

        }

        push(@peaks, $last);

        foreach my $peak (@peaks) {

            { lock(@qresults);
              push(@qresults, shared_clone([$id, $peak->{coords}->[0], $peak->{coords}->[1] + 1, sprintf("%." . $decimals . "f", mean(@{$peak->{score}})), sprintf("%." . $decimals . "e", pcombine($peak->{pvalue}, "S")), summit($peak->{coords}->[0], $peak->{coords}->[1], \@icov)])); }

            if ($metaplot) {

                my @metapeak = metacoverage([(0) x $peak->{coords}->[0], (1) x ($peak->{coords}->[1] - $peak->{coords}->[0] + 1), (0) x ($length - 1 - $peak->{coords}->[1])], SIZERNA);
                @pmeta = map { $pmeta[$_] + $metapeak[$_] } 0 .. SIZERNA - 1;

            }

            if ($metacoding) {

                my (@peakcov, @cov5utr, @cov3utr, @covcds);
                @peakcov = (0) x $length;
                $peakcov[$_] = 1 for ($peak->{coords}->[0] .. $peak->{coords}->[1]);
                @cov5utr = (0) x SIZE5UTR;
                @cov3utr = (0) x SIZE3UTR;

                @cov5utr = metacoverage([@peakcov[0 .. $orfStart - 1]], SIZE5UTR) if ($orfStart > 0);
                @covcds = metacoverage([@peakcov[$orfStart .. $orfEnd]], SIZECDS);
                @cov3utr = metacoverage([@peakcov[$orfEnd + 1 .. $length - 1]], SIZE3UTR) if ($orfEnd < $length - 1);

                undef(@peakcov);
                @peakcov = (@cov5utr, @covcds, @cov3utr);
                @pmetacod = map { $pmetacod[$_] + $peakcov[$_] } 0 .. $#peakcov;

                { lock(%results);
                  $results{codpeaks}++; }

            }

        }

        { lock(%results);
          $results{pass}++;
          $results{peaks} += @peaks; }

        if ($metaonlyenriched) {

            if ($control) {

                lock($ctrlnreads);
                $ctrlnreads += $cnreads;

            }

            { lock($ipnreads);
              $ipnreads += $ientry->readscount(); }

            if ($metaplot) {

                my (@metacoverage);

                if ($control) {

                    @metacoverage = metacoverage(\@ccov, SIZERNA);
                    @cmeta = map { $cmeta[$_] + $metacoverage[$_] } 0 .. SIZERNA - 1;

                    undef(@metacoverage);

                }

                @metacoverage = metacoverage(\@icov, SIZERNA);
                @imeta = map { $imeta[$_] + $metacoverage[$_] } 0 .. SIZERNA - 1;

            }

            if ($metacoding && $orfEnd) {

                my (@cmeta5, @cmeta3, @cmetac, @imeta5,
                    @imeta3, @imetac, @metacod);
                @cmeta5 = (0) x SIZE5UTR;
                @cmeta3 = (0) x SIZE3UTR;
                @imeta5 = @cmeta5;
                @imeta3 = @cmeta3;

                if ($orfStart > 0) {

                    @cmeta5 = metacoverage([@ccov[0 .. $orfStart - 1]], SIZE5UTR) if ($control);
                    @imeta5 = metacoverage([@icov[0 .. $orfStart - 1]], SIZE5UTR);

                }

                @cmetac = metacoverage([@ccov[$orfStart .. $orfEnd]], SIZECDS) if ($control);
                @imetac = metacoverage([@icov[$orfStart .. $orfEnd]], SIZECDS);

                if ($orfEnd < $length - 1) {

                    @cmeta3 = metacoverage([@ccov[$orfEnd + 1 .. $length - 1]], SIZE3UTR) if ($control);
                    @imeta3 = metacoverage([@icov[$orfEnd + 1 .. $length - 1]], SIZE3UTR);

                }

                if ($control) {

                    @metacod = (@cmeta5, @cmetac, @cmeta3);
                    $cmetacod[$_] += ($metacod[$_] ) for (0 .. $#metacod);

                    undef(@metacod);

                }

                @metacod = (@imeta5, @imetac, @imeta3);
                $imetacod[$_] += ($metacod[$_] ) for (0 .. $#metacod);

            }

        }

    }

    if ($metaplot) {

        if ($control) {

            lock(@metactrl);
            @metactrl = map { $metactrl[$_] + $cmeta[$_] } 0 .. SIZERNA - 1;

        }

        { lock(@metaip);
          @metaip = map { $metaip[$_] + $imeta[$_] } 0 .. SIZERNA - 1; }

        { lock(@metapeaks);
          @metapeaks = map { $metapeaks[$_] + $pmeta[$_] } 0 .. SIZERNA - 1; }

    }

    if ($metacoding) {

        if ($control) {

            lock(@metactrlcod);
            @metactrlcod = map { $metactrlcod[$_] + $cmetacod[$_] } 0 .. $codsize - 1;

        }

        { lock(@metaipcod);
          @metaipcod = map { $metaipcod[$_] + $imetacod[$_] } 0 .. $codsize - 1; }

        { lock(@metapeakscod);
          @metapeakscod = map { $metapeakscod[$_] + $pmetacod[$_] } 0 .. $codsize - 1; }

    }

}

sub cleanup {

    unlink(glob($tmpDir . "*"));

    rmtree($tmpDir);

}

sub refinepeak {

    my %parameters = @_;

    my ($start, $end, $imean, $cmean,
        $peakicovmedian);
    ($start, $end) = ($parameters{peak}->{coords}->[0], $parameters{peak}->{coords}->[1]);
    $peakicovmedian = median(@{$parameters{icov}});

    for (my $i = $parameters{peak}->{coords}->[0]; $i <= $parameters{peak}->{coords}->[1]; $i++) {

        my $enrichment = $relaxed ? mround($parameters{enrichment}->[$i], 0.5) : $parameters{enrichment}->[$i];

        if ($enrichment >= $enrich &&
            $parameters{icov}->[$i] >= $peakicovmedian) {

            $start = $i;

            last;

        }

    }

    if ($start == $parameters{peak}->{coords}->[0]) {

        for (my $i = $parameters{peak}->{coords}->[0]; $i >= 0; $i--) {

            my $enrichment = $relaxed ? mround($parameters{enrichment}->[$i], 0.5) : $parameters{enrichment}->[$i];

            last if ($enrichment < $enrich ||
                     $parameters{icov}->[$i] < $peakicovmedian);

            $start = $i;

        }

    }

    for (my $i = $parameters{peak}->{coords}->[1]; $i >= $parameters{peak}->{coords}->[0]; $i--) {

        my $enrichment = $relaxed ? mround($parameters{enrichment}->[$i], 0.5) : $parameters{enrichment}->[$i];

        if ($enrichment >= $enrich  &&
            $parameters{icov}->[$i] >= $peakicovmedian) {

            $end = $i;

            last;

        }

    }

    if ($end == $parameters{peak}->{coords}->[1]) {

        for (my $i = $parameters{peak}->{coords}->[1]; $i < $parameters{length}; $i++) {

            my $enrichment = $relaxed ? mround($parameters{enrichment}->[$i], 0.5) : $parameters{enrichment}->[$i];

            last if ($enrichment < $enrich ||
                     $parameters{icov}->[$i] < $peakicovmedian);

            $end = $i;

        }

    }

    ($parameters{peak}->{coords}->[0], $parameters{peak}->{coords}->[1]) = ($start, $end);

    $imean = mean(@{$parameters{icov}}[$start .. $end]);

    if ($control) {

        $cmean = mean(@{$parameters{ccov}}[$start .. $end]);
        @{$parameters{peak}->{score}} = logarithm((($imean + $pseudocount) / ($parameters{imedian} + $pseudocount)) / (($cmean + $pseudocount) / ($parameters{cmedian} + $pseudocount)), 2);

    }
    else {  @{$parameters{peak}->{score}} = logarithm((($imean + $pseudocount) / ($parameters{imedian} + $pseudocount)), 2); }

    return($parameters{peak});

}

sub summit {

    my ($start, $end, $icov) = @_;

    if ($summit &&
        $icov) {

        my ($max, $first, @max, @win, @last);
        $max = max(@{$icov}[$start .. $end]);
        @max = grep { $icov->[$_] == $max } $start .. $end;
        $first = shift(@max);
        @last = ($first, $first + 1);
        @win = @last;

        while(my $base = shift(@max)) {

            if ($base == $win[1]) { $win[1] = $base + 1; }
            else {

                if (($win[1] - $win[0]) > ($last[1] - $last[0])) {

                    @last = @win;
                    @win = ($base, $base + 1);

                }

            }

        }

        @last = @win if (($win[1] - $win[0]) > ($last[1] - $last[0]));

        return(@last);

    }

}

sub metacoverage {

    my ($coverage, $size) = @_;

    $size--;
    my @metacoverage = map { $coverage->[round($_ / $size * $#{$coverage})] } 0 .. $size;

    return(@metacoverage);

}

sub plotcoverage {

    my ($file, $ip, $ctrl, $peak) = @_;

    my ($maxy, $metaprofile, $metacodprofile);
    $metaprofile = $file =~ m/meta-gene/ ? 1 : 0;
    $metacodprofile = $file =~ m/meta-coding/ ? 1 : 0;
    $maxy = $control ? max(@{$ip}, @{$ctrl}) : max(@{$ip});
    $ctrl ||= [];
    $peak ||= [];

    if ($maxy) {

        my ($hasPeaks, $plot, $profilePlot, @plots, %legendColors, %xBreaks);
        $hasPeaks = ($metacodprofile || $metaprofile) && @$peak && sum(@$peak) ? 1 : 0;
        %legendColors = ("IP" => "#376eae");
        $legendColors{"Control"} = "#006400" if (@$ctrl);

        if ($metacodprofile) {

            %xBreaks = ( "0"    => "TSS",
                         "75"   => "Start",
                         "725"  => "Stop",
                         "1000" => "TES" );

        }
        elsif ($metaprofile) {

            %xBreaks = ( "0"    => "0",
                         "250"  => "25",
                         "500"  => "50",
                         "750"  => "75",
                         "1000" => "100" );

        }

        $plot = Graphics::Image->new( file   => $file,
                                      width  => 10,
                                      height => $hasPeaks ? 8 : 4,
                                      R      => $R,
                                      tmpdir => $tmpDir );

        $profilePlot = Graphics::Chart::Linedot->new( data           => [ @$ip, @$ctrl ],
                                                      dataLabels     => { sample => [ ("IP") x scalar(@$ip), ("Control") x scalar(@$ctrl) ],
                                                                          pos    => [ (1 .. scalar(@$ip)), (1 .. scalar(@$ctrl))] },
                                                      dataLabelType  => { pos => "numeric" },
                                                      x              => "pos",
                                                      fill           => "sample",
                                                      plotDataPoints => 0,
                                                      legendColors   => \%legendColors,
                                                      xTitle         => $metaprofile ? "Relative position (\%)" : ($metacodprofile ? undef : "Position (nt)"),
                                                      yTitle         => "RPM",
                                                      axisTitleSize  => 14,
                                                      labelTextSize  => 12,
                                                      legendKeyWidth => 14,
                                                      legendTextSize => 12,
                                                      xBreaks        => \%xBreaks );

        push(@plots, $profilePlot);

        if ($hasPeaks) {

            my $peakPlot = Graphics::Chart::Linedot->new( data           => $peak,
                                                          dataLabels     => { sample => [ ("Peaks") x scalar(@$peak) ],
                                                                              pos    => [ (1 .. scalar(@$peak)) ] },
                                                          dataLabelType  => { pos => "numeric" },
                                                          x              => "pos",
                                                          fill           => "sample",
                                                          plotDataPoints => 0,
                                                          legend         => 0,
                                                          legendColors   => { "Peaks" => "#000000" },
                                                          xTitle         => $metaprofile ? "Relative position (\%)" : undef,
                                                          yTitle         => "% Peaks",
                                                          axisTitleSize  => 14,
                                                          legendKeyWidth => 14,
                                                          legendTextSize => 12,
                                                          labelTextSize  => 12,
                                                          xBreaks        => \%xBreaks );

            push(@plots, $peakPlot);

        }

        $plot->plot([ map { [ $_ ] } @plots ]);

        # if ($control) {

        #     my $ctrlPlot = Graphics::Chart::Linedot->new( x              => "fp",
        #                                                   legend         => 1,
        #                                                   data           => \@tp,
        #                                                   dataLabels     => { fp => \@fp,
        #                                                                       id => \@id },
        #                                                   dataLabelType  => { fp => "numeric"},
        #                                                   legendSort     => [ (grep { $_ ne "Overall" } sort(uniq(@id))), "Overall" ],
        #                                                   fill           => "id",
        #                                                   xLabelAngle    => 90,
        #                                   plotDataPoints => 0,
        #                                   colorPalette   => "Paired",
        #                                   xTitle         => "False Positive Rate",
        #                                   yTitle         => "True Positive Rate" );
        #     $ccov = Graphics::Object::Path->new( fontsize    => 15,
        #                                           height      => 1,
        #                                           yrange      => [0, $maxy],
        #                                           forceyrange => 1,
        #                                           values      => $ctrl,
        #                                           stroke      => "rgb(0, 100, 0)",
        #                                           yname       => "Control (RPM)" );

        # }

        # $icov = Graphics::Object::Path->new( fontsize    => 15,
        #                                      height      => 1,
        #                                      yrange      => [0, $maxy],
        #                                      forceyrange => 1,
        #                                      values      => $ip,
        #                                      stroke      => "rgb(55, 110, 174)",
        #                                      yname       => "IP (RPM)" );

        # $ruler = Graphics::Object::Ruler->new( fontsize => 15,
        #                                        height   => 1,
        #                                        labelMap => $metacodprofile ? { "0" => "TSS",
        #                                                                        "7.5" => "Start",
        #                                                                        "72.5" => "Stop",
        #                                                                        "100"  => "TES" } : {},
        #                                        range    => $metaprofile || $metacodprofile ? [0, 100] : [0, scalar(@{$ip})],
        #                                        ticksby  => $metaprofile ? 10 : max(1, round(scalar(@{$ip}) / 10)),
        #                                        name     => $metaprofile ? "Relative position (\%)" : ($metacodprofile ? undef : "Position (nt)"),
        #                                        toend    => 1 );

        # $image->addobjects($ccov) if ($control);
        # $image->addobject($icov);

        # if (($metaprofile || $metacodprofile) &&
        #     @{$peak}) {

        #     $pcov = Graphics::Object::Path->new( fontsize    => 15,
        #                                          height      => 1,
        #                                          values      => $peak,
        #                                          stroke      => "rgb(0, 0, 0)",
        #                                          yname       => "\% of peaks" );

        #     $image->addobject($pcov);

        # }

        # $image->addobject($ruler);

        # $io->write($image->xml());

    }

}

sub help {

    print "\n  [!] Error: Invalid option. Please check the help\n" if ($_[0]);

    die <<HELP;

 RF PeakCall (v$Core::Utils::VERSION)
 RNA Framework [http://www.rnaframework.com]

 Author:  Danny Incarnato (dincarnato[at]rnaframework.com)
 Summary: Performs peak calling of RNA immunoprecipitation (RIP) experiments

 Usage:   rf-peakcall [Options]

 Options                                            Description
 -c   or --control              <string>             Control sample RC file
 -I   or --IP                   <string>             Immunoprecipitated (IP) sample RC file
 -i   or --index                <string>[,<string>]  A comma separated (no spaces) list of RCI index files for the provided RC files.
                                                     Note: RCI files must be provided in the order: 1. Control, 2. IP.
                                                     If a single RCI file is specified along with both control and IP sample,
                                                     it will be used for all samples.
                                                     If no RCI index is provided, it will be created at runtime, and stored in
                                                     the same folder of the control, and IP samples.
 -l   or --whitelist            <string>             A whitelist containing transcript IDs (one per each row) to restrict the analysis to
 -p   or --processors           <int>                Number of processors to use (Default: 1)
 -o   or --output               <string>             Output folder (Default: <IP>_vs_<Control> if a control RC file is provided, or
                                                                             <IP> if only the IP RC file is provided)
 -ow  or --overwrite                                 Overwrites output folder (if the specified folder already exists)

 Peak calling options
 -w   or --window               <int>                Window's size for peak calling (>=10, Default: 150)
 -f   or --offset               <int>                Offset for window sliding (>=1, Default: Window's size / 2)
 -md  or --merge-distance       <int>                Maximum distance in nt for merging non-overlapping windows (>=0, Default: 50)
 -e   or --enrichment           <float>              Minimum log2 enrichment in IP vs. control sample, for reporting a peak (>=1, Default: 3)
 -v   or --p-value              <float>              P-value cutoff for reporting a peak (0<=p<=1, Default: 0.05)
 -s   or --summit                                    Generates an additional BED file containing the coordinates of peak summits
 -r   or --refine                                    Refines peak boundaries
 -x   or --relaxed                                   Uses more relaxed criteria to refine peak boundaries (requires -r)
 -pc  or --pseudocount          <float>              Pseudocount added to read counts to avoid division by 0 (>0, Default: 1)
 -mc  or --mean-coverage        <float>              Discards any transcript with mean coverage below this threshold in the control sample (>=0, Default: 0)
 -ec  or --median-coverage      <float>              Discards any transcript with median coverage below this threshold in the control sample (>=0, Default: 0)
 -D   or --decimals             <int>                Number of decimals for reporting enrichment/p-value (1-10, Default: 3)

 Plotting options 
 -g   or --img                                       Enables the generation of coverage plots (one per transcript; requires R)
 -mp  or --meta-plot                                 Enables the generation of meta-gene plots (requires R)
 -mcp or --meta-coding-plot                          Enables the generation of a protein-coding-only meta-gene plot, by aligning the TSS, start codon,
                                                     stop codon, and TES (requires R)
 -of  or --orf-file             <string>             Path to a BED file containing the transcript-level coordinates of the CDSs
                                                     Note: if no file is provided, the longest ORF will be automatically identified using the following parameters
 -mo  or --min-orf-length       <int>                Minimum length (in aa) to select the longest ORF (requires -mcp, Default: 50)
 -als or --alt-start                                 The longest ORF is allowed to start with alternative start codons (requires -mcp)
 -ans or --any-start                                 The longest ORF is allowed to start with any codon (requires -mcp)
 -gc  or --genetic-code         <int>                Genetic code table for the reference organism (requires -mcp, 1-33, Default: 1)
 -eo  or --plot-enriched-only                        Meta-gene plots will be calculated only on transcripts with a detected enrichment (requires -mp or -mcp)
 -R   or --R-path               <string>             Path to R executable (Default: assumes R is in PATH)

HELP

}
